BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKSLKK.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 1998 Caltech, UCSB
10//
11// Module: EvtD0ToKSLKK.cc
12//
13// Description: Routine to handle D0 -> KS0 K+ K- and KL0 K+ K-
14//
15// Modification history:
16//
17// Liaoyuan Dong Oct 24 15:17:18 2023 Module created
18//
19//------------------------------------------------------------------------
20//
21#include "EvtD0ToKSLKK.hh"
22#include "../EvtGenBase/EvtConst.hh"
23#include "../EvtGenBase/EvtDecayTable.hh"
24#include "../EvtGenBase/EvtGenKine.hh"
25#include "../EvtGenBase/EvtPDL.hh"
26#include "../EvtGenBase/EvtParticle.hh"
27#include "../EvtGenBase/EvtPatches.hh"
28#include "../EvtGenBase/EvtReport.hh"
29#include "../EvtGenBase/EvtResonance.hh"
30#include "../EvtGenBase/EvtResonance2.hh"
31#include <fstream>
32#include <stdlib.h>
33#include <string>
34using std::endl;
35
36EvtD0ToKSLKK::~EvtD0ToKSLKK() {}
37
38void EvtD0ToKSLKK::getName( std::string& model_name ) { model_name = "D0ToKSLKK"; }
39
40EvtDecayBase* EvtD0ToKSLKK::clone() { return new EvtD0ToKSLKK; }
41
42void EvtD0ToKSLKK::init() {
43 // checkNArg(1);
44 Narg = getNArg();
45 if ( Narg == 0 ) { Uspin = 0; }
46 else { Uspin = getArg( 0 ); }
47 charge = EvtPDL::getStdHep( getParentId() );
48 Daug0Id = EvtPDL::getStdHep( getDaug( 0 ) );
49 checkNDaug( 3 );
50 checkSpinParent( EvtSpinType::SCALAR );
51 checkSpinDaughter( 0, EvtSpinType::SCALAR );
52 checkSpinDaughter( 1, EvtSpinType::SCALAR );
53 checkSpinDaughter( 2, EvtSpinType::SCALAR );
54
55 phi[0] = 0;
56 rho[0] = 1;
57 phi[1] = 3.367900972;
58 rho[1] = 0.5567227875;
59 phi[2] = 3.568151995;
60 rho[2] = 1.122904371;
61 phi[3] = 3.353326851;
62 rho[3] = 1.414616255;
63 phi[4] = -0.7616302559;
64 rho[4] = 1.623197367;
65
66 modetype[0] = 23;
67 modetype[1] = 23;
68 modetype[2] = 13;
69 modetype[3] = 13;
70 modetype[4] = 13;
71
72 if ( Uspin == 1 && Daug0Id == 130 )
73 {
74 modetype[0] = 232;
75 modetype[1] = 232;
76 modetype[2] = 13;
77 modetype[3] = 13;
78 modetype[4] = 13;
79 }
80
81 /*
82 for (int i=0; i<5; i++) {
83 cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
84 }
85 */
86 width[0] = 0.272;
87 width[1] = 0.004249;
88 width[2] = 0.272;
89 width[3] = 0.258;
90 width[4] = 0.4;
91
92 mass[0] = 0.919;
93 mass[1] = 1.019461;
94 mass[2] = 0.919;
95 mass[3] = 1.439;
96 mass[4] = 1.439;
97
98 mDM = 1.86484;
99 mK0 = 0.497614;
100 mKa = 0.49368;
101 mPi = 0.13957;
102 mEta = 0.547862;
103 mKa2 = 0.24371994; // 0.49368^2;
104 mPi2 = 0.01947978; // 0.13957^2;
105 mEta2 = 0.30015277; // 0.547862^2;
106 mass_EtaP = 0.95778;
107 mass_Kaon = 0.49368;
108 mass_KS = 0.4976;
109
110 math_pi = 3.1415926;
111 mass_Pion2 = 0.0194797849;
112 mass_2Pion = 0.27914;
113 math_2pi = 6.2831852;
114 rD2 = 25.0; // 5*5
115 rRes2 = 9.0; // 3*3
116 g2 = 0.23; // K*0(1430)
117
118 GS1 = 0.636619783;
119 GS2 = 0.01860182466;
120 GS3 = 0.1591549458; // 1/(2*math_2pi)
121 GS4 = 0.00620060822; // mass_Pion2/math_pi
122
123 rho_omega = 0.00294;
124 phi_omega = -0.02;
125
126 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
127 for ( int i = 0; i < 4; i++ )
128 {
129 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
130 }
131}
132
133void EvtD0ToKSLKK::initProbMax() {
134 setProbMax( 1700.0 );
135 if ( Daug0Id == 310 ) setProbMax( 1700.0 );
136 else if ( Daug0Id == 130 && Uspin == 1 ) setProbMax( 1785.0 );
137 else setProbMax( 1700.0 );
138}
139
140void EvtD0ToKSLKK::decay( EvtParticle* p ) {
141 /*
142 double maxprob = 0.0;
143 for(int ir=0;ir<=60000000;ir++){
144
145 p->initializePhaseSpace(getNDaug(),getDaugs());
146 EvtVector4R D1 = p->getDaug(0)->getP4();
147 EvtVector4R D2 = p->getDaug(1)->getP4();
148 EvtVector4R D3 = p->getDaug(2)->getP4();
149
150 //charge = EvtPDL::getStdHep(p->getId());
151
152 double P1[4], P2[4], P3[4];
153 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
154 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
155 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
156
157 if(Daug0Id==310) SorL = true; else SorL = false;
158 double value;
159 int spin[5]={0,1,0,0,1};
160 if(SorL){
161 int g0[5]={5,1,3,1,5};
162 double r0[5] = {3,3,3,3,3};
163 double r1[5] = {5,5,5,5,5};
164 int nstates=5;
165 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, r0, r1,
166 value, 0, nstates,charge,SorL); }else if((!SorL)&&Uspin==1){ int g0[5]={5,1,3,1,5}; double
167 r0[5] = {-1.566394443,-1.33043736,3,3,3}; double r1[5] =
168 {0.1844175671,-1.397710917,5,5,5}; int nstates=5; calEva(P1, P2, P3, mass, width, rho,
169 phi, g0, spin, modetype, r0, r1, value, 0, nstates,charge,SorL); }else{ int
170 g0[5]={5,1,3,1,5}; double r0[5] = {3,3,3,3,3}; double r1[5] = {5,5,5,5,5}; int nstates=5;
171 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, r0, r1,
172 value, 0, nstates,charge,SorL);
173 }
174
175 if (value<0) continue;
176 if(value>maxprob) {
177 maxprob=value;
178 cout << "ir= " << ir << endl;
179 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3]
180 <<"};"<< endl; cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<<
181 P2[3] <<"};"<< endl; cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2]
182 <<","<< P3[3] <<"};"<< endl; cout << "MAX====> " << maxprob << endl;
183 }
184 }
185 printf("MAXprob = %.10f\n",maxprob);
186 */
187
188 p->initializePhaseSpace( getNDaug(), getDaugs() );
189 EvtVector4R D1 = p->getDaug( 0 )->getP4();
190 EvtVector4R D2 = p->getDaug( 1 )->getP4();
191 EvtVector4R D3 = p->getDaug( 2 )->getP4();
192
193 // charge = EvtPDL::getStdHep(p->getId());
194
195 double P1[4], P2[4], P3[4];
196 P1[0] = D1.get( 0 );
197 P1[1] = D1.get( 1 );
198 P1[2] = D1.get( 2 );
199 P1[3] = D1.get( 3 );
200 P2[0] = D2.get( 0 );
201 P2[1] = D2.get( 1 );
202 P2[2] = D2.get( 2 );
203 P2[3] = D2.get( 3 );
204 P3[0] = D3.get( 0 );
205 P3[1] = D3.get( 1 );
206 P3[2] = D3.get( 2 );
207 P3[3] = D3.get( 3 );
208
209 // P1[0] = 0.662259; P1[1] = -0.277017; P1[2] = 0.211748; P1[3] = -0.263429;
210 // P2[0] = 0.664422; P2[1] = 0.283792; P2[2] = -0.243945; P2[3] = 0.240194;
211 // P3[0] = 0.561591; P3[1] = 0.188743; P3[2] = -0.0572267;P3[3] = -0.181021;
212
213 // P1[0] = 0.663036; P1[1] = -0.188126; P1[2] = 0.265033; P1[3] = -0.293883;
214 // P2[0] = 0.694463; P2[1] = 0.201378; P2[2] = -0.399611; P2[3] = -0.195755;
215 // P3[0] = 0.530119; P3[1] = 0.122859; P3[2] = -0.123733; P3[3] = -0.083098;
216
217 if ( Daug0Id == 310 ) SorL = true;
218 else SorL = false;
219 double value;
220 int spin[5] = { 0, 1, 0, 0, 1 };
221 if ( SorL )
222 {
223 int g0[5] = { 5, 1, 3, 1, 5 };
224 double r0[5] = { 3, 3, 3, 3, 3 };
225 double r1[5] = { 5, 5, 5, 5, 5 };
226 int nstates = 5;
227 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, r0, r1, value, 0, nstates,
228 charge, SorL );
229 }
230 else if ( ( !SorL ) && Uspin == 1 )
231 {
232 int g0[5] = { 5, 1, 3, 1, 5 };
233 double r0[5] = { -1.566394443, -1.33043736, 3, 3, 3 };
234 double r1[5] = { 0.1844175671, -1.397710917, 5, 5, 5 };
235 int nstates = 5;
236 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, r0, r1, value, 0, nstates,
237 charge, SorL );
238 }
239 else
240 {
241 int g0[5] = { 5, 1, 3, 1, 5 };
242 double r0[5] = { 3, 3, 3, 3, 3 };
243 double r1[5] = { 5, 5, 5, 5, 5 };
244 int nstates = 5;
245 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, r0, r1, value, 0, nstates,
246 charge, SorL );
247 }
248
249 setProb( value );
250
251 return;
252}
253
254void EvtD0ToKSLKK::Com_Multi( double a1[2], double a2[2], double res[2] ) {
255 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
256 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
257}
258void EvtD0ToKSLKK::Com_Divide( double a1[2], double a2[2], double res[2] ) {
259 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
260 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
261 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
262}
263//------------base---------------------------------
264double EvtD0ToKSLKK::SCADot( double a1[4], double a2[4] ) {
265 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
266 return _cal;
267}
268double EvtD0ToKSLKK::barrier( int l, double sa, double sb, double sc, double r, double mass ) {
269 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
270 // if(q < 0) q = 1e-16;
271 if ( q < 0 ) q = -q;
272 double z;
273 z = q * r * r;
274 double sa0;
275 sa0 = mass * mass;
276 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
277 // if(q0 < 0) q0 = 1e-16;
278 if ( q0 < 0 ) q0 = -q0;
279 double z0 = q0 * r * r;
280 double F = 0.0;
281 if ( l == 0 ) F = 1;
282 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
283 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
284 return F;
285}
286//------------------spin-------------------------------------------
287void EvtD0ToKSLKK::calt1( double daug1[4], double daug2[4], double t1[4] ) {
288 double p, pq, tmp;
289 double pa[4], qa[4];
290 for ( int i = 0; i < 4; i++ )
291 {
292 pa[i] = daug1[i] + daug2[i];
293 qa[i] = daug1[i] - daug2[i];
294 }
295 p = SCADot( pa, pa );
296 pq = SCADot( pa, qa );
297 tmp = pq / p;
298 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
299}
300void EvtD0ToKSLKK::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
301 double p, r;
302 double pa[4], t1[4];
303 calt1( daug1, daug2, t1 );
304 r = SCADot( t1, t1 ) / 3.0;
305 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
306 p = SCADot( pa, pa );
307 for ( int i = 0; i < 4; i++ )
308 {
309 for ( int j = 0; j < 4; j++ )
310 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
311 }
312}
313//-------------------prop--------------------------------------------
314void EvtD0ToKSLKK::propagatorCBW( double mass, double width, double sx, double prop[2] ) {
315 double a[2], b[2];
316 a[0] = 1;
317 a[1] = 0;
318 b[0] = mass * mass - sx;
319 b[1] = -mass * width;
320 Com_Divide( a, b, prop );
321}
322double EvtD0ToKSLKK::wid( double mass2, double mass, double sa, double sb, double sc,
323 double r2, int l ) {
324 double widm = 0.;
325 double m = sqrt( sa );
326 double tmp = sb - sc;
327 double tmp1 = sa + tmp;
328 double q = 0.25 * tmp1 * tmp1 / sa - sb;
329 // if(q<0) q = 1e-16;
330 if ( q < 0 ) q = -q;
331 double tmp2 = mass2 + tmp;
332 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
333 // if(q0<0) q0 = 1e-16;
334 if ( q0 < 0 ) q0 = -q0;
335 double z = q * r2;
336 double z0 = q0 * r2;
337 double t = q / q0;
338 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
339 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
340 else if ( l == 2 )
341 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
342 return widm;
343}
344double EvtD0ToKSLKK::widl1( double mass2, double mass, double sa, double sb, double sc,
345 double r2 ) {
346 double widm = 0.;
347 double m = sqrt( sa );
348 double tmp = sb - sc;
349 double tmp1 = sa + tmp;
350 double q = 0.25 * tmp1 * tmp1 / sa - sb;
351 // if(q<0) q = 1e-16;
352 if ( q < 0 ) q = -q;
353 double tmp2 = mass2 + tmp;
354 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
355 // if(q0<0) q0 = 1e-16;
356 if ( q0 < 0 ) q0 = -q0;
357 double z = q * r2;
358 double z0 = q0 * r2;
359 double F = ( 1 + z0 ) / ( 1 + z );
360 double t = q / q0;
361 widm = t * sqrt( t ) * mass / m * F;
362 return widm;
363}
364void EvtD0ToKSLKK::propagatorRBW( double mass, double width, double sa, double sb, double sc,
365 double r2, int l, double prop[2] ) {
366 double a[2], b[2];
367 double mass2 = mass * mass;
368
369 a[0] = 1;
370 a[1] = 0;
371 b[0] = mass2 - sa;
372 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
373 Com_Divide( a, b, prop );
374}
375
376void EvtD0ToKSLKK::propagatorFlatte( double mass, double width, double sa, double prop[2] ) {
377
378 double q2_Pi, q2_Ka;
379 double rhoPi[2], rhoKa[2];
380
381 q2_Pi = 0.25 * sa - mPi * mPi;
382 q2_Ka = 0.25 * sa - mKa * mKa;
383
384 if ( q2_Pi > 0 )
385 {
386 rhoPi[0] = 2.0 * sqrt( q2_Pi / sa );
387 rhoPi[1] = 0.0;
388 }
389 if ( q2_Pi <= 0 )
390 {
391 rhoPi[0] = 0.0;
392 rhoPi[1] = 2.0 * sqrt( -q2_Pi / sa );
393 }
394
395 if ( q2_Ka > 0 )
396 {
397 rhoKa[0] = 2.0 * sqrt( q2_Ka / sa );
398 rhoKa[1] = 0.0;
399 }
400 if ( q2_Ka <= 0 )
401 {
402 rhoKa[0] = 0.0;
403 rhoKa[1] = 2.0 * sqrt( -q2_Ka / sa );
404 }
405
406 // from paper PLB 598(2004) 149 : The simga pole in J/psi -> omega pi+ pi-
407 // Double_t g1 = 0.138;
408 // Double_t g2 = 0.6141;// g1*4.45;
409 // M = 0.970
410
411 // PLB 607(2005) 243 : Resonances in J/psi -> phi pi+ pi- and phi K+ K-
412 // M = 0.965
413 // Double_t g1 = 0.165;
414 // Double_t g2 = 0.69465;// g1*4.21;
415
416 // from paper PRD 86 052006 (2012) :LHCb barB^0_s -> J/psi pi+pi-
417 // Double_t g1 = 0.199;
418 // Double_t g2 = 0.5970;// g1*3.0;
419
420 double a[2], b[2];
421 a[0] = 1;
422 a[1] = 0;
423 b[0] = mass * mass - sa + 0.165 * rhoPi[1] + 0.69465 * rhoKa[1];
424 b[1] = -( 0.165 * rhoPi[0] + 0.69465 * rhoKa[0] );
425 Com_Divide( a, b, prop );
426}
427
428//------------GS---used by rho----------------------------
429void EvtD0ToKSLKK::propagatorGS( double mass, double width, double sa, double sb, double sc,
430 double r2, double prop[2] ) {
431 double a[2], b[2];
432 double mass2 = mass * mass;
433 double tmp = sb - sc;
434 double tmp1 = sa + tmp;
435 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
436 // if(q2<0) q2 = 1e-16;
437 if ( q2 < 0 ) q2 = -q2;
438
439 double tmp2 = mass2 + tmp;
440 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
441 // if(q02<0) q02 = 1e-16;
442 if ( q02 < 0 ) q02 = -q02;
443
444 double q = sqrt( q2 );
445 double q0 = sqrt( q02 );
446 double m = sqrt( sa );
447 double q03 = q0 * q02;
448 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
449
450 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
451 double h0 = GS1 * q0 / mass * tmp3;
452 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
453 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
454 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
455
456 a[0] = 1.0 + d * width / mass;
457 a[1] = 0.0;
458 b[0] = mass2 - sa + width * f;
459 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
460 Com_Divide( a, b, prop );
461}
462void EvtD0ToKSLKK::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
463 const double m1430 = 1.441;
464 const double sa0 = 2.076481; // m1430*m1430;
465 const double w1430 = 0.193;
466 const double Lass1 = 0.25 / sa0;
467 double tmp = sb - sc;
468 double tmp1 = sa0 + tmp;
469 double q0 = Lass1 * tmp1 * tmp1 - sb;
470 // if(q0<0) q0 = 1e-16;
471 if ( q0 < 0 ) q0 = -q0;
472 double tmp2 = sa + tmp;
473 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
474 double q = sqrt( qs );
475 double width = w1430 * q * m1430 / sqrt( sa * q0 );
476 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
477 if ( temp_R < 0 ) temp_R += math_pi;
478 double deltaR = -109.7 * math_pi / 180.0 + temp_R;
479 double temp_F =
480 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
481 if ( temp_F < 0 ) temp_F += math_pi;
482 double deltaF = 0.1 * math_pi / 180.0 + temp_F;
483 double deltaS = deltaR + 2.0 * deltaF;
484 double t1 = 0.96 * sin( deltaF );
485 double t2 = sin( deltaR );
486 double CF[2], CS[2];
487 CF[0] = cos( deltaF );
488 CF[1] = sin( deltaF );
489 CS[0] = cos( deltaS );
490 CS[1] = sin( deltaS );
491 prop[0] = t1 * CF[0] + t2 * CS[0];
492 prop[1] = t1 * CF[1] + t2 * CS[1];
493}
494
495void EvtD0ToKSLKK::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
496 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
497 if ( q > 0 )
498 {
499 rho[0] = 2 * sqrt( q / sa );
500 rho[1] = 0;
501 }
502 else if ( q < 0 )
503 {
504 rho[0] = 0;
505 rho[1] = 2 * sqrt( -q / sa );
506 }
507}
508
509void EvtD0ToKSLKK::propagatorKstr1430( double mass, double sx, double* sb, double* sc,
510 double prop[2] ) // K*1430 Flatte
511{
512 double unit[2] = { 1.0 };
513 double ci[2] = { 0, 1 };
514 double rho1[2];
515 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
516 double rho2[2];
517 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
518 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
519 double tmp1[2] = { gKPi_Kstr1430, 0 };
520 double tmp11[2];
521 double tmp2[2] = { gEtaPK_Kstr1430, 0 };
522 double tmp22[2];
523 Com_Multi( tmp1, rho1, tmp11 );
524 Com_Multi( tmp2, rho2, tmp22 );
525 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
526 double tmp31[2];
527 Com_Multi( tmp3, ci, tmp31 );
528 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
529 Com_Divide( unit, tmp4, prop );
530}
531
532void EvtD0ToKSLKK::propagatora0980p( double mass, double sx, double* sb, double* sc,
533 double prop[2] ) // a0980p Flatte
534{
535 double unit[2] = { 1.0 };
536 double ci[2] = { 0, 1 };
537 double rho1[2];
538 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
539 double rho2[2];
540 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
541 double gKK_a0980p = 0.341 * 0.892, gEtapi_a0980p = 0.341;
542 double tmp1[2] = { gKK_a0980p, 0 };
543 double tmp11[2];
544 double tmp2[2] = { gEtapi_a0980p, 0 };
545 double tmp22[2];
546 Com_Multi( tmp1, rho1, tmp11 );
547 Com_Multi( tmp2, rho2, tmp22 );
548 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
549 double tmp31[2];
550 Com_Multi( tmp3, ci, tmp31 );
551 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
552 Com_Divide( unit, tmp4, prop );
553}
554
555void EvtD0ToKSLKK::propagatora0980pfloated( double mass, double sx, double* sb, double* sc,
556 double gKK, double prop[2] ) // a0980p Flatte
557{
558 double unit[2] = { 1.0 };
559 double ci[2] = { 0, 1 };
560 double rho1[2];
561 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
562 double rho2[2];
563 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
564 double gKK_a0980p = 0.341 * gKK, gEtapi_a0980p = 0.341;
565 double tmp1[2] = { gKK_a0980p, 0 };
566 double tmp11[2];
567 double tmp2[2] = { gEtapi_a0980p, 0 };
568 double tmp22[2];
569 Com_Multi( tmp1, rho1, tmp11 );
570 Com_Multi( tmp2, rho2, tmp22 );
571 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
572 double tmp31[2];
573 Com_Multi( tmp3, ci, tmp31 );
574 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
575 Com_Divide( unit, tmp4, prop );
576}
577
578void EvtD0ToKSLKK::propagatora0980wm( double mass, double width, double sx, double sb,
579 double sc, double prop[2] ) // a0980wm Flatte
580{
581 double unit[2] = { 1.0 };
582 double ci[2] = { 0, 1 };
583 double rho1[2];
584 double tmp1[2] = { mass, 0 };
585 double tmp2[2] = { width, 0 };
586 double tmp11[2], tmp22[2];
587 Flatte_rhoab( sx, sb, sc, rho1 );
588 Com_Multi( tmp1, rho1, tmp11 );
589 Com_Multi( tmp11, tmp2, tmp22 );
590 double tmp3[2];
591 Com_Multi( tmp22, ci, tmp3 );
592 double tmp4[2] = { mass * mass - sx - tmp3[0], -1.0 * tmp3[1] };
593 Com_Divide( unit, tmp4, prop );
594}
595
596void EvtD0ToKSLKK::propagatora09800( double mass, double sx, double* sb, double* sc,
597 double prop[2] ) // a09800 Flatte
598{
599 double unit[2] = { 1.0 };
600 double ci[2] = { 0, 1 };
601 double rho1[2];
602 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
603 double rho2[2];
604 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
605 double rho3[2];
606 Flatte_rhoab( sx, sb[2], sc[2], rho3 );
607 double gKK_a09800 = 0.341 * 0.892, gEtapi_a09800 = 0.341, gK0K0_a09800 = 0.341 * 0.892;
608 double tmp1[2] = { gKK_a09800, 0 };
609 double tmp11[2];
610 double tmp2[2] = { gEtapi_a09800, 0 };
611 double tmp22[2];
612 double tmp3[2] = { gK0K0_a09800, 0 };
613 double tmp33[2];
614 Com_Multi( tmp1, rho1, tmp11 );
615 Com_Multi( tmp2, rho2, tmp22 );
616 Com_Multi( tmp3, rho3, tmp33 );
617 double tmp4[2] = { tmp11[0] + tmp22[0] + tmp33[0], tmp11[1] + tmp22[1] + tmp33[1] };
618 double tmp41[2];
619 Com_Multi( tmp4, ci, tmp41 );
620 double tmp5[2] = { mass * mass - sx - tmp41[0], -1.0 * tmp41[1] };
621 Com_Divide( unit, tmp5, prop );
622}
623
624void EvtD0ToKSLKK::propagatora09800floated( double mass, double sx, double* sb, double* sc,
625 double gKK, double prop[2] ) // a09800 Flatte
626{
627 double unit[2] = { 1.0 };
628 double ci[2] = { 0, 1 };
629 double rho1[2];
630 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
631 double rho2[2];
632 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
633 double rho3[2];
634 Flatte_rhoab( sx, sb[2], sc[2], rho3 );
635 double gKK_a09800 = 0.341 * gKK, gEtapi_a09800 = 0.341, gK0K0_a09800 = 0.341 * gKK;
636 double tmp1[2] = { gKK_a09800, 0 };
637 double tmp11[2];
638 double tmp2[2] = { gEtapi_a09800, 0 };
639 double tmp22[2];
640 double tmp3[2] = { gK0K0_a09800, 0 };
641 double tmp33[2];
642 Com_Multi( tmp1, rho1, tmp11 );
643 Com_Multi( tmp2, rho2, tmp22 );
644 Com_Multi( tmp3, rho3, tmp33 );
645 double tmp4[2] = { tmp11[0] + tmp22[0] + tmp33[0], tmp11[1] + tmp22[1] + tmp33[1] };
646 double tmp41[2];
647 Com_Multi( tmp4, ci, tmp41 );
648 double tmp5[2] = { mass * mass - sx - tmp41[0], -1.0 * tmp41[1] };
649 Com_Divide( unit, tmp5, prop );
650}
651
652void EvtD0ToKSLKK::propagatora098002channel( double mass, double sx, double* sb, double* sc,
653 double prop[2] ) // a09800 Flatte
654{
655 double unit[2] = { 1.0 };
656 double ci[2] = { 0, 1 };
657 double rho1[2];
658 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
659 double rho2[2];
660 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
661 double gKK_a09800 = 0.341 * 0.892, gEtapi_a09800 = 0.341;
662 double tmp1[2] = { gKK_a09800, 0 };
663 double tmp11[2];
664 double tmp2[2] = { gEtapi_a09800, 0 };
665 double tmp22[2];
666 Com_Multi( tmp1, rho1, tmp11 );
667 Com_Multi( tmp2, rho2, tmp22 );
668 double tmp4[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
669 double tmp41[2];
670 Com_Multi( tmp4, ci, tmp41 );
671 double tmp5[2] = { mass * mass - sx - tmp41[0], -1.0 * tmp41[1] };
672 Com_Divide( unit, tmp5, prop );
673}
674
675void EvtD0ToKSLKK::rhoab( double sa, double sb, double sc, double res[2] ) {
676 double tmp = sa + sb - sc;
677 double q = 0.25 * tmp * tmp / sa - sb;
678 if ( q >= 0 )
679 {
680 res[0] = 2.0 * sqrt( q / sa );
681 res[1] = 0.0;
682 }
683 else
684 {
685 res[0] = 0.0;
686 res[1] = 2.0 * sqrt( -q / sa );
687 }
688}
689void EvtD0ToKSLKK::rho4Pi( double sa, double res[2] ) {
690 double temp = 1.0 - 0.3116765584 / sa; // 0.3116765584=0.13957*0.13957*16
691 if ( temp >= 0 )
692 {
693 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
694 res[1] = 0.0;
695 }
696 else
697 {
698 res[0] = 0.0;
699 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) ); //?????????????????????
700 }
701}
702
703void EvtD0ToKSLKK::propagatorsigma500( double sa, double sb, double sc,
704 double prop[2] ) { // for pipi_s wave
705 double f = 0.5843 + 1.6663 * sa;
706 const double M = 0.9264;
707 const double mass2 = 0.85821696; // M*M
708 const double mpi2d2 = 0.00973989245;
709 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
710 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
711 rhoab( sa, sb, sc, rho1s );
712 rhoab( mass2, sb, sc, rho1M );
713 rho4Pi( sa, rho2s );
714 rho4Pi( mass2, rho2M );
715 Com_Divide( rho1s, rho1M, rho1 );
716 Com_Divide( rho2s, rho2M, rho2 );
717 double a[2], b[2];
718 a[0] = 1.0;
719 a[1] = 0.0;
720 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
721 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
722 Com_Divide( a, b, prop );
723}
724
725void EvtD0ToKSLKK::getprop( double sa, double sb, double sc, double mass, double width,
726 double prop[2] ) { // rho propagator
727 double prop1[2], prop2[2];
728 propagatorGS( mass, width, sa, sb, sc, 9.0, prop1 );
729 propagatorRBW( 0.783, 0.008, sa, sb, sc, 3.0, 1, prop2 );
730 double coef_omega[2];
731 coef_omega[0] = rho_omega * cos( phi_omega ), coef_omega[1] = rho_omega * sin( phi_omega );
732 double one[2];
733 one[0] = 1;
734 one[1] = 0;
735 double temp[2];
736 Com_Multi( coef_omega, prop2, temp );
737 temp[0] = one[0] + 0.783 * 0.783 * temp[0];
738 temp[1] = one[1] + 0.783 * 0.783 * temp[1];
739 Com_Multi( prop1, temp, prop );
740}
741double EvtD0ToKSLKK::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
742 double mass ) {
743 double pR[4], pD[4];
744 double temp_PDF, v_re;
745 temp_PDF = 0.0;
746 v_re = 0.0;
747 double B[2], s1, s2, s3, sR, sD;
748 for ( int i = 0; i < 4; i++ )
749 {
750 pR[i] = P1[i] + P2[i];
751 pD[i] = pR[i] + P3[i];
752 }
753 s1 = SCADot( P1, P1 );
754 s2 = SCADot( P2, P2 );
755 s3 = SCADot( P3, P3 );
756 sR = SCADot( pR, pR );
757 sD = SCADot( pD, pD );
758 int G[4][4];
759 for ( int i = 0; i != 4; i++ )
760 {
761 for ( int j = 0; j != 4; j++ )
762 {
763 if ( i == j )
764 {
765 if ( i == 0 ) G[i][j] = 1;
766 else G[i][j] = -1;
767 }
768 else G[i][j] = 0;
769 }
770 }
771 if ( Ang == 0 )
772 {
773 B[0] = 1;
774 B[1] = 1;
775 temp_PDF = 1;
776 }
777 if ( Ang == 1 )
778 {
779 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
780 B[1] = barrier( 1, sD, sR, s3, 5.0, mDM );
781 // B[0] = Barrier(1,sR,s1,s2,9.0);
782 // B[1] = Barrier(1,sD,sR,s3,25.0);
783 double t1[4], T1[4];
784 calt1( P1, P2, t1 );
785 calt1( pR, P3, T1 );
786 temp_PDF = 0;
787 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
788 }
789 if ( Ang == 2 )
790 {
791 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
792 B[1] = barrier( 2, sD, sR, s3, 5.0, mDM );
793 // B[0] = Barrier(2,sR,s1,s2,9.0);
794 // B[1] = Barrier(2,sD,sR,s3,25.0);
795 double t2[4][4], T2[4][4];
796 calt2( P1, P2, t2 );
797 calt2( pR, P3, T2 );
798 temp_PDF = 0;
799 for ( int i = 0; i < 4; i++ )
800 {
801 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
802 }
803 }
804 v_re = temp_PDF * B[0] * B[1];
805 return v_re;
806}
807
808void EvtD0ToKSLKK::calPDF( double* Ks0, double* K1, double* K2, double* mass1, double* width1,
809 double* amp, double* phase, int* g0, int* spin, int* modetype,
810 double* r0, double* r1, int first, int last, double PDF[2] ) {
811 double P12[4], P23[4], P13[4];
812 double cof[2], amp_PDF[2];
813 double s12, s13, s23;
814 for ( int i = 0; i < 4; i++ )
815 {
816 P12[i] = K1[i] + Ks0[i];
817 P13[i] = K2[i] + Ks0[i];
818 P23[i] = K1[i] + K2[i];
819
820 // printf("abcde: %lf, %lf, %lf, %lf, %lf\n",P12[i], P13[i], K1[i], K2[i], Ks0[i]);
821 }
822 s12 = SCADot( P12, P12 );
823 s13 = SCADot( P13, P13 );
824 s23 = SCADot( P23, P23 );
825 double s1, s2, s3;
826 s1 = SCADot( Ks0, Ks0 );
827 s2 = SCADot( K1, K1 );
828 s3 = SCADot( K2, K2 );
829 double pro[2], temp_PDF, amp_tmp[2];
830 double Amp_KPiS[2];
831 // double mass1sq;
832 amp_PDF[0] = 0;
833 amp_PDF[1] = 0;
834 PDF[0] = 0;
835 PDF[1] = 0;
836 amp_tmp[0] = 0;
837 amp_tmp[1] = 0;
838
839 // printf("ssssss: %lf, %lf, %lf, %lf, %lf, %lf\n",s1,s2,s3,s12,s13,s23);
840
841 for ( int i = first; i < last; i++ )
842 {
843 amp_tmp[0] = 0;
844 amp_tmp[1] = 0;
845 // mass1sq = mass1[i]*mass1[i];
846 cof[0] = amp[i] * cos( phase[i] );
847 cof[1] = amp[i] * sin( phase[i] );
848 temp_PDF = 0;
849
850 if ( modetype[i] == 12 )
851 {
852 temp_PDF = DDalitz( Ks0, K1, K2, spin[i], mass1[i] );
853 if ( g0[i] == 1 )
854 propagatorRBW( mass1[i], width1[i], s12, mass_KS * mass_KS, mKa2, rRes2, spin[i],
855 pro );
856 if ( g0[i] == 2 )
857 { // a0980p Flatte
858 double s11[2] = { mass_KS * mass_KS, mPi * mPi };
859 double s22[2] = { mKa * mKa, mEta * mEta };
860 propagatora0980p( mass1[i], s12, s11, s22, pro );
861 }
862 if ( g0[i] == 3 )
863 propagatora0980wm( mass1[i], width1[i], s12, mass_KS * mass_KS, mKa * mKa, pro );
864 if ( g0[i] == 0 )
865 {
866 pro[0] = 1;
867 pro[1] = 0;
868 }
869 amp_tmp[0] = temp_PDF * pro[0];
870 amp_tmp[1] = temp_PDF * pro[1];
871 }
872
873 if ( modetype[i] == 13 )
874 {
875 temp_PDF = DDalitz( Ks0, K2, K1, spin[i], mass1[i] );
876 if ( g0[i] == 1 )
877 propagatorRBW( mass1[i], width1[i], s13, mass_KS * mass_KS, mKa2, rRes2, spin[i],
878 pro );
879 if ( g0[i] == 2 )
880 { // a0980p Flatte
881 double s11[2] = { mass_KS * mass_KS, mPi * mPi };
882 double s22[2] = { mKa * mKa, mEta * mEta };
883 propagatora0980p( mass1[i], s13, s11, s22, pro );
884 }
885 if ( g0[i] == 3 )
886 propagatora0980wm( mass1[i], width1[i], s13, mass_KS * mass_KS, mKa * mKa, pro );
887 if ( g0[i] == 4 )
888 { // a0980p Flatte
889 double s11[2] = { mass_KS * mass_KS, mPi * mPi };
890 double s22[2] = { mKa * mKa, mEta * mEta };
891 propagatora0980pfloated( mass1[i], s13, s11, s22, r0[i], pro );
892 }
893 if ( g0[i] == 5 )
894 propagatorGS( mass1[i], width1[i], s13, mass_KS * mass_KS, mKa2, rRes2, pro );
895 if ( g0[i] == 0 )
896 {
897 pro[0] = 1;
898 pro[1] = 0;
899 }
900 // printf("%lf, %lf\n",pro[0],pro[1]);
901 amp_tmp[0] = temp_PDF * pro[0];
902 amp_tmp[1] = temp_PDF * pro[1];
903 }
904
905 if ( modetype[i] == 23 )
906 {
907 temp_PDF = DDalitz( K1, K2, Ks0, spin[i], mass1[i] );
908 if ( g0[i] == 0 )
909 {
910 pro[0] = 1;
911 pro[1] = 0;
912 }
913 if ( g0[i] == 1 )
914 propagatorRBW( mass1[i], width1[i], s23, mKa2, mKa2, rRes2, spin[i], pro );
915 if ( g0[i] == 2 )
916 { // a09800 Flatte
917 double s11[3] = { mKa * mKa, mPi * mPi, mass_KS * mass_KS };
918 double s22[3] = { mKa * mKa, mEta * mEta, mass_KS * mass_KS };
919 propagatora09800( mass1[i], s23, s11, s22, pro );
920 }
921 if ( g0[i] == 3 )
922 { // a09800 Flatte 2 channel
923 double s11[3] = { mKa * mKa, mPi * mPi };
924 double s22[3] = { mKa * mKa, mEta * mEta };
925 propagatora098002channel( mass1[i], s23, s11, s22, pro );
926 }
927 if ( g0[i] == 4 )
928 { propagatorFlatte( mass1[i], width1[i], s23, pro ); } // Only for f0(980)
929 if ( g0[i] == 5 )
930 propagatora0980wm( mass1[i], width1[i], s23, mKa * mKa, mKa * mKa, pro );
931 if ( g0[i] == 6 )
932 { // a09800 Flatte
933 double s11[3] = { mKa * mKa, mPi * mPi, mass_KS * mass_KS };
934 double s22[3] = { mKa * mKa, mEta * mEta, mass_KS * mass_KS };
935 propagatora09800floated( mass1[i], s23, s11, s22, r0[i], pro );
936 }
937 amp_tmp[0] = temp_PDF * pro[0];
938 amp_tmp[1] = temp_PDF * pro[1];
939 }
940
941 //------------------------------------KLKK-------------------------------------------------------
942 if ( modetype[i] == 232 )
943 {
944 temp_PDF = DDalitz( K1, K2, Ks0, spin[i], mass1[i] );
945 if ( g0[i] == 0 )
946 {
947 pro[0] = 1;
948 pro[1] = 0;
949 }
950 if ( g0[i] == 1 )
951 propagatorRBW( mass1[i], width1[i], s23, mKa2, mKa2, rRes2, spin[i], pro );
952 if ( g0[i] == 2 )
953 { // a09800 Flatte
954 double s11[3] = { mKa * mKa, mPi * mPi, mass_KS * mass_KS };
955 double s22[3] = { mKa * mKa, mEta * mEta, mass_KS * mass_KS };
956 propagatora09800( mass1[i], s23, s11, s22, pro );
957 }
958 if ( g0[i] == 3 )
959 { // a09800 Flatte 2 channel
960 double s11[3] = { mKa * mKa, mPi * mPi };
961 double s22[3] = { mKa * mKa, mEta * mEta };
962 propagatora098002channel( mass1[i], s23, s11, s22, pro );
963 }
964 if ( g0[i] == 4 )
965 { propagatorFlatte( mass1[i], width1[i], s23, pro ); } // Only for f0(980)
966 if ( g0[i] == 5 )
967 propagatora0980wm( mass1[i], width1[i], s23, mKa * mKa, mKa * mKa, pro );
968 if ( g0[i] == 6 )
969 { // a09800 Flatte
970 double s11[3] = { mKa * mKa, mPi * mPi, mass_KS * mass_KS };
971 double s22[3] = { mKa * mKa, mEta * mEta, mass_KS * mass_KS };
972 propagatora09800floated( mass1[i], s23, s11, s22, r0[i], pro );
973 }
974 if ( g0[i] == 7 ) propagatorGS( mass1[i], width1[i], s23, mKa * mKa, mKa2, rRes2, pro );
975
976 double uspinbr[2], coff[2];
977 coff[0] = r0[i] * cos( r1[i] );
978 coff[1] = r0[i] * sin( r1[i] );
979 double unit[2];
980 unit[0] = 1.0;
981 unit[1] = 0.0;
982
983 uspinbr[0] =
984 unit[0] - 2 * ( ( 0.22650 * 0.22650 ) / ( 1. - ( 0.22650 * 0.22650 ) ) ) * coff[0];
985 uspinbr[1] =
986 unit[1] - 2 * ( ( 0.22650 * 0.22650 ) / ( 1. - ( 0.22650 * 0.22650 ) ) ) * coff[1];
987
988 double amp_tmpp[2];
989 amp_tmpp[0] = temp_PDF * pro[0];
990 amp_tmpp[1] = temp_PDF * pro[1];
991
992 Com_Multi( uspinbr, amp_tmpp, amp_tmp );
993 }
994
995 if ( modetype[i] == 100 )
996 {
997 amp_tmp[0] = 1.0;
998 amp_tmp[1] = 0.0;
999 }
1000
1001 // if(modetype[i] == 132){
1002 // KPiSLASS(s13,s1,s3,Amp_KPiS);
1003 // amp_tmp[0] = Amp_KPiS[0];
1004 // amp_tmp[1] = Amp_KPiS[1];
1005 // }
1006
1007 Com_Multi( amp_tmp, cof, amp_PDF );
1008 PDF[0] += amp_PDF[0];
1009 PDF[1] += amp_PDF[1];
1010 // printf("%.10lf, %.10lf\n",PDF[0],PDF[1]);
1011 }
1012}
1013
1014void EvtD0ToKSLKK::calEva( double* Ks0, double* K1, double* K2, double* mass1, double* width1,
1015 double* amp, double* phase, int* g0, int* spin, int* modetype,
1016 double* r0, double* r1, double& Result, int first, int last,
1017 int charge, bool SorL ) {
1018 double PDFD0[2], PDFD0bar[2];
1019 if ( charge > 0 )
1020 {
1021 calPDF( Ks0, K1, K2, mass1, width1, amp, phase, g0, spin, modetype, r0, r1, first, last,
1022 PDFD0 );
1023 Ks0[0] = Ks0[0];
1024 Ks0[1] = ( -1.0 ) * Ks0[1];
1025 Ks0[2] = ( -1.0 ) * Ks0[2];
1026 Ks0[3] = ( -1.0 ) * Ks0[3];
1027 K1[0] = K1[0];
1028 K1[1] = ( -1.0 ) * K1[1];
1029 K1[2] = ( -1.0 ) * K1[2];
1030 K1[3] = ( -1.0 ) * K1[3];
1031 K2[0] = K2[0];
1032 K2[1] = ( -1.0 ) * K2[1];
1033 K2[2] = ( -1.0 ) * K2[2];
1034 K2[3] = ( -1.0 ) * K2[3];
1035
1036 // calPDF(Ks0,K2,K1,mass1,width1,amp,phase,g0,spin,modetype,r0,r1,first,last,PDFD0bar);
1037 calPDF( Ks0, K1, K2, mass1, width1, amp, phase, g0, spin, modetype, r0, r1, first, last,
1038 PDFD0bar );
1039 }
1040 if ( charge < 0 )
1041 {
1042 calPDF( Ks0, K1, K2, mass1, width1, amp, phase, g0, spin, modetype, r0, r1, first, last,
1043 PDFD0bar );
1044 Ks0[0] = Ks0[0];
1045 Ks0[1] = ( -1.0 ) * Ks0[1];
1046 Ks0[2] = ( -1.0 ) * Ks0[2];
1047 Ks0[3] = ( -1.0 ) * Ks0[3];
1048 K1[0] = K1[0];
1049 K1[1] = ( -1.0 ) * K1[1];
1050 K1[2] = ( -1.0 ) * K1[2];
1051 K1[3] = ( -1.0 ) * K1[3];
1052 K2[0] = K2[0];
1053 K2[1] = ( -1.0 ) * K2[1];
1054 K2[2] = ( -1.0 ) * K2[2];
1055 K2[3] = ( -1.0 ) * K2[3];
1056 // calPDF(Ks0,K2,K1,mass1,width1,amp,phase,g0,spin,modetype,r0,r1,first,last,PDFD0);
1057 calPDF( Ks0, K1, K2, mass1, width1, amp, phase, g0, spin, modetype, r0, r1, first, last,
1058 PDFD0 );
1059 }
1060 /*
1061 //-----------------------Quantum Correlation
1062 Correction--------------------------------- double r_tag = 0.0586; double R_tag = 1;
1063 double delta_tag = 192.1/180.0*3.1415926;
1064 double qcf[2];
1065 qcf[0] = r_tag*R_tag*cos(-1.0*delta_tag);
1066 qcf[1] = r_tag*R_tag*sin(-1.0*delta_tag);
1067 double ampD0_part1[2], qcfampD0bar[2];
1068 Com_Multi(qcf,PDFD0bar,qcfampD0bar);
1069 if(SorL){
1070 ampD0_part1[0] = PDFD0[0] - qcfampD0bar[0];
1071 ampD0_part1[1] = PDFD0[1] - qcfampD0bar[1];
1072 } else{
1073 ampD0_part1[0] = PDFD0[0] + qcfampD0bar[0];
1074 ampD0_part1[1] = PDFD0[1] + qcfampD0bar[1];
1075 }
1076 double value = ampD0_part1[0]*ampD0_part1[0]+ampD0_part1[1]*ampD0_part1[1] +
1077 r_tag*r_tag*(1-R_tag*R_tag)*(PDFD0bar[0]*PDFD0bar[0]+PDFD0bar[1]*PDFD0bar[1]);
1078 */
1079 double value = PDFD0[0] * PDFD0[0] + PDFD0[1] * PDFD0[1];
1080 if ( value <= 0 ) value = 1e-20;
1081 Result = value;
1082
1083 // printf("value: %lf\n",value);
1084}
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
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
one
const DifNumber one
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
EvtD0ToKSLKK::clone
EvtDecayBase * clone()
Definition EvtD0ToKSLKK.cc:40
EvtD0ToKSLKK::~EvtD0ToKSLKK
virtual ~EvtD0ToKSLKK()
Definition EvtD0ToKSLKK.cc:36
EvtD0ToKSLKK::getName
void getName(std::string &name)
Definition EvtD0ToKSLKK.cc:38
EvtD0ToKSLKK::decay
void decay(EvtParticle *p)
Definition EvtD0ToKSLKK.cc:140
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::getArg
double getArg(int j)
Definition EvtDecayBase.cc:520
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:261
EvtDecayBase::getParentId
EvtId getParentId()
Definition EvtDecayBase.hh:63
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:465
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:68
EvtDecayBase::getDaug
EvtId getDaug(int i)
Definition EvtDecayBase.hh:69
EvtDecayProb::setProb
void setProb(double prob)
Definition EvtDecayProb.hh:40
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition EvtPDL.hh:61
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