BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSKSpi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of the EvtGen package developed jointly
4// for the BaBar and CLEO collaborations. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/COPYRIGHT
8// Copyright (C) 1998 Caltech, UCSB
9//
10// Module: EvtDsToKSKSpi.cc
11//
12// Description: Routine to handle three-body decays of Ds+/Ds-
13//
14// Modification history:
15// Liaoyuan Dong Sun Jun 02 01:23:25 2024 Module created
16//------------------------------------------------------------------------
18#include <fstream>
19#include <stdlib.h>
22#include "EvtGenBase/EvtPDL.hh"
27#include <string>
30using std::endl;
31
33
34void EvtDsToKSKSpi::getName(std::string& model_name){
35 model_name="DsToKSKSpi";
36}
37
41
43 // check that there are 0 arguments
44 checkNArg(0);
45 checkNDaug(3);
50
51 phi[0] = 0;
52 rho[0] = 1;
53 phi[1] = 0;
54 rho[1] = 0;
55 phi[2] = 0;
56 rho[2] = 0;
57 phi[3] = 0;
58 rho[3] = 0;
59 phi[4] = 0;
60 rho[4] = 0;
61 phi[0] = 0.0;
62 phi[1] = 2.3351;
63
64 rho[0] = 1.0;
65 rho[1] = 3.3891;
66
67 modetype[0]= 13;
68 modetype[1]= 12;
69 modetype[2]= 12;
70 modetype[3]= 13;
71 modetype[4]= 13;
72
73 //cout << "DsToKSKSpi :" << endl;
74 //for (int i=0; i<5; i++) {
75 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
76 //}
77
78 width[0] = 5.0800e-02;
79 width[1] = 1.4041e-01;//1.2300e-01
80 width[2] = 2.6500e-01;
81 width[3] = 2.7000e-01;
82 width[4] = 5.0800e-02;
83
84 mass[0] = 8.9166e-01;
85 mass[1] = 1.7228e+00;
86 mass[2] = 1.3500e+00;
87 mass[3] = 1.4250e+00;
88 mass[4] = 8.9166e-01;
89
90 mD = 1.86486;
91 mDs = 1.9683;
92 rRes = 9.0;
93 rD = 5.0;
94 metap = 0.95778;
95 mkstr = 0.89594;
96 mk0 = 0.497614;
97 mass_Kaon = 0.49368;
98 mass_Pion = 0.13957;
99 mass_Pi0 = 0.1349766;
100 math_pi = 3.1415926;
101 //mrho = 0.77549;
102 //Grho = 0.1491;
103 ma0 = 0.99;
104 Ga0 = 0.0756;
105 meta = 0.547862;
106
107 GS1 = 0.636619783;
108 GS2 = 0.01860182466;
109 GS3 = 0.1591549458; // 1/(2*math_2pi)
110 GS4 = 0.00620060822; // mass_Pion2/math_pi
111
112 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
113 for (int i=0; i<4; i++) {
114 for (int j=0; j<4; j++) {
115 G[i][j] = GG[i][j];
116 }
117 }
118}
119
121 setProbMax(1156.0);
122}
123
125/*
126 double maxprob = 0.0;
127 for(int ir=0;ir<=60000000;ir++){
128 p->initializePhaseSpace(getNDaug(),getDaugs());
129 EvtVector4R D1 = p->getDaug(0)->getP4();
130 EvtVector4R D2 = p->getDaug(1)->getP4();
131 EvtVector4R D3 = p->getDaug(2)->getP4();
132
133 double P1[4], P2[4], P3[4];
134 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
135 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
136 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
137
138 double value;
139// value = calDalEva(P1, P2, P3);
140 int g0[5]={1,1,1,1,12};
141 int spin[5]={1,0,0,0,0};
142 int nstates=5;
143 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
144 if (value<0) continue;
145 if(value>maxprob) {
146 maxprob=value;
147 cout << "ir= " << ir << endl;
148 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
149 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
150 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
151 cout << "MAX====> " << maxprob << endl;
152 }
153 }
154 printf("MAXprob = %.10f\n",maxprob);
155*/
156
158 EvtVector4R D1 = p->getDaug(0)->getP4();
159 EvtVector4R D2 = p->getDaug(1)->getP4();
160 EvtVector4R D3 = p->getDaug(2)->getP4();
161
162 double P1[4], P2[4], P3[4];
163 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
164 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
165 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
166// P1[0] = 0.770687885731853; P1[1] = 0.085867043334332; P1[2] =-0.550183164280213; P1[3] =-0.190434925445587;
167// P2[0] = 0.840768850984605; P2[1] =-0.240025188001374; P2[2] = 0.631078772181167; P2[3] = 0.058310035302289;
168// P3[0] = 0.393776666168180; P3[1] =-0.159532664043002; P3[2] =-0.300922058149019; P3[3] = 0.139912371490140;
169
170 double value;
171 int g0[5]={1,1,1,1,12};
172 int spin[5]={1,0,0,0,0};
173 int nstates=5;
174 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
175
176 setProb(value);
177 return ;
178
179}
180
181
182
183void EvtDsToKSKSpi::Com_Multi(double a1[2], double a2[2], double res[2])
184{
185 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
186 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
187}
188void EvtDsToKSKSpi::Com_Divide(double a1[2], double a2[2], double res[2])
189{
190 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
191 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
192 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
193}
194//------------base---------------------------------
195double EvtDsToKSKSpi::SCADot(double a1[4], double a2[4])
196{
197 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
198 return _cal;
199}
200double EvtDsToKSKSpi::barrier(int l, double sa, double sb, double sc, double r, double mass)
201{
202 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
203 if(q < 0) q = 1e-16;
204 double z;
205 z = q*r*r;
206 double sa0;
207 sa0 = mass*mass;
208 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
209 if(q0 < 0) q0 = 1e-16;
210 double z0 = q0*r*r;
211 double F = 0.0;
212 if(l == 0) F = 1;
213 if(l == 1) F = sqrt((1+z0)/(1+z));
214 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
215 return F;
216}
217
218void EvtDsToKSKSpi::calt1(double daug1[4], double daug2[4], double t1[4])
219{
220 double p, pq, tmp;
221 double pa[4], qa[4];
222 for(int i=0; i<4; i++) {
223 pa[i] = daug1[i] + daug2[i];
224 qa[i] = daug1[i] - daug2[i];
225 }
226 p = SCADot(pa,pa);
227 pq = SCADot(pa,qa);
228 tmp = pq/p;
229 for(int i=0; i<4; i++) {
230 t1[i] = qa[i] - tmp*pa[i];
231 }
232}
233void EvtDsToKSKSpi::calt2(double daug1[4], double daug2[4], double t2[4][4])
234{
235 double p, r;
236 double pa[4], t1[4];
237 calt1(daug1,daug2,t1);
238 r = SCADot(t1,t1)/3.0;
239 for(int i=0; i<4; i++) {
240 pa[i] = daug1[i] + daug2[i];
241 }
242 p = SCADot(pa,pa);
243 for(int i=0; i<4; i++) {
244 for(int j=0; j<4; j++) {
245 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
246 }
247 }
248}
249//-------------------prop--------------------------------------------
250
251double EvtDsToKSKSpi::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
252{
253 double widm = 0.;
254 double m = sqrt(sa);
255 double tmp = sb-sc;
256 double tmp1 = sa+tmp;
257 double q = 0.25*tmp1*tmp1/sa-sb;
258 if(q<0) q = 1e-16;
259 double tmp2 = mass2+tmp;
260 double q0 = 0.25*tmp2*tmp2/mass2-sb;
261 if(q0<0) q0 = 1e-16;
262 double z = q*r2;
263 double z0 = q0*r2;
264 double t = q/q0;
265 if(l == 0) {widm = sqrt(t)*mass/m;}
266 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
267 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
268 return widm;
269}
270double EvtDsToKSKSpi::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
271{
272 double widm = 0.;
273 double m = sqrt(sa);
274 double tmp = sb-sc;
275 double tmp1 = sa+tmp;
276 double q = 0.25*tmp1*tmp1/sa-sb;
277 if(q<0) q = 1e-16;
278 double tmp2 = mass2+tmp;
279 double q0 = 0.25*tmp2*tmp2/mass2-sb;
280 if(q0<0) q0 = 1e-16;
281 double z = q*r2;
282 double z0 = q0*r2;
283 double F = (1+z0)/(1+z);
284 double t = q/q0;
285 widm = t*sqrt(t)*mass/m*F;
286 return widm;
287}
288void EvtDsToKSKSpi::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
289{
290 double a[2], b[2];
291 a[0] = 1;
292 a[1] = 0;
293 b[0] = mass2-sa;
294 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
295 Com_Divide(a,b,prop);
296}
297
298void EvtDsToKSKSpi::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
299{
300 double a[2], b[2];
301 double tmp = sb-sc;
302 double tmp1 = sa+tmp;
303 double q2 = 0.25*tmp1*tmp1/sa-sb;
304 if(q2<0) q2 = 1e-16;
305
306 double tmp2 = mass2+tmp;
307 double q02 = 0.25*tmp2*tmp2/mass2-sb;
308 if(q02<0) q02 = 1e-16;
309
310 double q = sqrt(q2);
311 double q0 = sqrt(q02);
312 double m = sqrt(sa);
313 double q03 = q0*q02;
314 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
315
316 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
317 double h0 = GS1*q0/mass*tmp3;
318 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
319 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
320 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
321
322 a[0] = 1.0+d*width/mass;
323 a[1] = 0.0;
324 b[0] = mass2-sa+width*f;
325 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
326 Com_Divide(a,b,prop);
327}
328
329
330void EvtDsToKSKSpi::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
331 const double m1430 = 1.441;
332 const double sa0 = 1.441*1.441; // m1430*m1430;
333 const double w1430 = 0.193;
334 const double Lass1 = 0.25/sa0;
335 double tmp = sb-sc;
336 double tmp1 = sa0+tmp;
337 double q0 = Lass1*tmp1*tmp1-sb;
338 if(q0<0) q0 = 1e-16;
339 double tmp2 = sa+tmp;
340 double qs = 0.25*tmp2*tmp2/sa-sb;
341 double q = sqrt(qs);
342 double width = w1430*q*m1430/sqrt(sa*q0);
343 double temp_R = atan(m1430*width/(sa0-sa));
344 if(temp_R<0) temp_R += math_pi;
345 double deltaR = -109.7 + temp_R;
346 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
347 if(temp_F<0) temp_F += math_pi;
348 double deltaF = 0.1 + temp_F;
349 double deltaS = deltaR + 2.0*deltaF;
350 double t1 = 0.96*sin(deltaF);
351 double t2 = sin(deltaR);
352 double CF[2], CS[2];
353 CF[0] = cos(deltaF);
354 CF[1] = sin(deltaF);
355 CS[0] = cos(deltaS);
356 CS[1] = sin(deltaS);
357 prop[0] = t1*CF[0] + t2*CS[0];
358 prop[1] = t1*CF[1] + t2*CS[1];
359}
360
361double EvtDsToKSKSpi::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
362 double pR[4], pD[4];
363 double temp_PDF, v_re;
364 temp_PDF = 0.0;
365 v_re = 0.0;
366 double B[2], s1, s2, s3, sR, sD;
367 for(int i=0; i<4; i++){
368 pR[i] = P1[i] + P2[i];
369 pD[i] = pR[i] + P3[i];
370 }
371 s1 = SCADot(P1,P1);
372 s2 = SCADot(P2,P2);
373 s3 = SCADot(P3,P3);
374 sR = SCADot(pR,pR);
375 sD = SCADot(pD,pD);
376 //int G[4][4];
377 //for(int i=0; i!=4; i++){
378 // for(int j=0; j!=4; j++){
379 // if(i==j){
380 // if(i==0) G[i][j] = 1;
381 // else G[i][j] = -1;
382 // }
383 // else G[i][j] = 0;
384 // }
385 //}
386 if(Ang == 0){
387 B[0] = 1;
388 B[1] = 1;
389 temp_PDF = 1;
390 }
391 if(Ang == 1){
392 B[0] = barrier(1,sR,s1,s2,3.0,mass);
393 B[1] = barrier(1,sD,sR,s3,5.0,1.9683);
394 double t1[4], T1[4];
395 calt1(P1,P2,t1);
396 calt1(pR,P3,T1);
397 temp_PDF = 0;
398 for(int i=0; i<4; i++){
399 temp_PDF += t1[i]*T1[i]*G[i][i];
400 }
401 }
402 if(Ang == 2){
403 B[0] = barrier(2,sR,s1,s2,3.0,mass);
404 B[1] = barrier(2,sD,sR,s3,5.0,1.9683);
405 double t2[4][4], T2[4][4];
406 calt2(P1,P2,t2);
407 calt2(pR,P3,T2);
408 temp_PDF = 0;
409 for(int i=0; i<4; i++){
410 for(int j=0; j<4; j++){
411 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
412 }
413 }
414 }
415 v_re = temp_PDF*B[0]*B[1];
416 return v_re;
417}
418
419
420void EvtDsToKSKSpi::calEva(double* Ks01, double* Ks02, double* Pic, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result)
421{
422 //double Ks0[4], Pic[4], Pi0[4];
423 //double rRes = 9.0;
424 double P12[4], P23[4], P13[4];
425 double cof[2], amp_PDF[2], PDF[2];
426 double scpi, sks02, sks01;
427 double s12, s13, s23;
428 for(int i=0; i<4; i++){
429 P12[i] = Ks01[i] + Ks02[i];
430 P13[i] = Ks01[i] + Pic[i];
431 P23[i] = Ks02[i] + Pic[i];
432 }
433 scpi = SCADot(Pic,Pic);
434 sks01 = SCADot(Ks01,Ks01);
435 sks02 = SCADot(Ks02,Ks02);
436 s12 = SCADot(P12,P12);
437 s13 = SCADot(P13,P13);
438 s23 = SCADot(P23,P23);
439 double pro[2], temp_PDF, amp_tmp[2],temp_PDF1 ,temp_PDF2,pro1[2],pro2[2];
440 double mass1sq;
441 amp_PDF[0] = 0;
442 amp_PDF[1] = 0;
443 PDF[0] = 0;
444 PDF[1] = 0;
445 amp_tmp[0] = 0;
446 amp_tmp[1] = 0;
447 for(int i=0; i<nstates; i++) {
448 amp_tmp[0] = 0;
449 amp_tmp[1] = 0;
450 mass1sq = mass1[i]*mass1[i];
451 cof[0] = amp[i]*cos(phase[i]);
452 cof[1] = amp[i]*sin(phase[i]);
453 temp_PDF = 0;
454/* if(modetype[i] == 23){
455 temp_PDF = DDalitz( Ks02, Pi0, Ks0, 1, mass1[i]);
456 if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],s23,scpi,snpi,9.0,pro);
457 if(g0[i]==0){
458 pro[0] = 1;
459 pro[1] = 0;
460 }
461 amp_tmp[0] = temp_PDF*pro[0];
462 amp_tmp[1] = temp_PDF*pro[1];
463
464 }
465*/
466 if(modetype[i] == 12){
467 temp_PDF = DDalitz(Ks01, Ks02, Pic, spin[i], mass1[i]);
468 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,sks01,sks02,rRes,spin[i],pro);
469 if(g0[i]==12) KPiSLASS(s12,sks01,sks02,pro);
470 if(g0[i]==0){
471 pro[0] = 1;
472 pro[1] = 0;
473 }
474 amp_tmp[0] = temp_PDF*pro[0];
475 amp_tmp[1] = temp_PDF*pro[1];
476 }
477 if(modetype[i] == 13){
478 temp_PDF1 = DDalitz(Ks01, Pic, Ks02, spin[i], mass1[i]);
479 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,sks01,scpi,rRes,spin[i],pro1);
480 if(g0[i]==12) KPiSLASS(s13,sks01,scpi,pro1);
481 if(g0[i]==0){
482 pro1[0] = 1;
483 pro1[1] = 0;
484 }
485 temp_PDF2 = DDalitz(Ks02, Pic, Ks01, spin[i], mass1[i]);
486 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s23,sks02,scpi,rRes,spin[i],pro2);
487 if(g0[i]==12) KPiSLASS(s23,sks02,scpi,pro2);
488 if(g0[i]==0){
489 pro2[0] = 1;
490 pro2[1] = 0;
491}
492 amp_tmp[0] = ( temp_PDF1*pro1[0]+temp_PDF2*pro2[0]);
493 amp_tmp[1] = ( temp_PDF1*pro1[1]+temp_PDF2*pro2[1]);
494 }
495 Com_Multi(amp_tmp,cof,amp_PDF);
496 PDF[0] += amp_PDF[0];
497 PDF[1] += amp_PDF[1];
498 }
499 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
500 //cout<<"value = "<<value<<endl;
501 if(value <=0) value = 1e-20;
502 Result = value;
503 // cout<<"value = "<<value<<" Result = "<<Result<<endl;
504}
double sin(const BesAngle a)
Definition BesAngle.h:210
double cos(const BesAngle a)
Definition BesAngle.h:213
double mass
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
****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
TCrossPart * CS
Definition Mcgpj.cxx:51
TTree * t
Definition binning.cxx:23
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
void checkSpinParent(EvtSpinType::spintype sp)
void setProbMax(double prbmx)
void checkNDaug(int d1, int d2=-1)
EvtId * getDaugs()
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
void setProb(double prob)
void decay(EvtParticle *p)
virtual ~EvtDsToKSKSpi()
EvtDecayBase * clone()
void getName(std::string &name)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
double get(int i) const
double double double double * s12
Definition qcdloop1.h:77
double double double double double * s23
Definition qcdloop1.h:77
const double b
Definition slope.cxx:9