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