BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSpipi0.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: EvtDsToKSpipi0.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 EvtDsToKSpipi0::getName(std::string& model_name){
35 model_name="DsToKSpipi0";
36}
37
41
43 // check that there are 0 arguments
44 checkNArg(0);
45 checkNDaug(3);
50
51 phi[0] = 1.9740e-01;
52 phi[1] = 3.1956e+00;
53 phi[3] = 1.9523e-01;
54 phi[4] = 2.1602e+00;
55
56 rho[0] = 3.0648e-01;
57 rho[1] = 4.0708e-01;
58 rho[3] = 7.5484e-01;
59 rho[4] = 2.6511e+00;
60
61 phi[2] = 0;
62 rho[2] = 1;
63 modetype[0]= 12;
64 modetype[1]= 13;
65 modetype[2]= 23;
66 modetype[3]= 131;
67 modetype[4]= 231;
68
69 //cout << "DsToKSpipi0 :" << endl;
70 //for (int i=0; i<5; i++) {
71 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
72 //}
73
74 width[0] = 5.0800e-02;
75 width[1] = 4.7300e-02;
76 width[2] = 1.5020e-01;
77 width[3] = 2.3200e-01;
78 width[4] = 4.0000e-01;
79 mass[0] = 8.9166e-01;
80 mass[1] = 8.9555e-01;
81 mass[2] = 7.6650e-01;
82 mass[3] = 1.4140e+00;
83 mass[4] = 1.4650e+00;
84
85 mD = 1.86486;
86 mDs = 1.9683;
87 rRes = 9.0;
88 rD = 5.0;
89 metap = 0.95778;
90 mkstr = 0.89594;
91 mk0 = 0.497614;
92 mass_Kaon = 0.49368;
93 mass_Pion = 0.13957;
94 mass_Pi0 = 0.1349766;
95 math_pi = 3.1415926;
96 //mrho = 0.77549;
97 //Grho = 0.1491;
98 ma0 = 0.99;
99 Ga0 = 0.0756;
100 meta = 0.547862;
101
102 GS1 = 0.636619783;
103 GS2 = 0.01860182466;
104 GS3 = 0.1591549458; // 1/(2*math_2pi)
105 GS4 = 0.00620060822; // mass_Pion2/math_pi
106
107 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
108 for (int i=0; i<4; i++) {
109 for (int j=0; j<4; j++) {
110 G[i][j] = GG[i][j];
111 }
112 }
113}
114
116 setProbMax(310.0);
117}
118
120 /*
121 double maxprob = 0.0;
122 for(int ir=0;ir<=6000000;ir++){
123 p->initializePhaseSpace(getNDaug(),getDaugs());
124 EvtVector4R D1 = p->getDaug(0)->getP4();
125 EvtVector4R D2 = p->getDaug(1)->getP4();
126 EvtVector4R D3 = p->getDaug(2)->getP4();
127
128 double P1[4], P2[4], P3[4];
129 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
130 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
131 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
132
133 double value;
134 //value = calDalEva(P1, P2, P3);
135 int g0[5]={1,1,1,1,1};
136 int nstates=5;
137 calEva(P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value);
138 if (value<0) continue;
139 if(value>maxprob) {
140 maxprob=value;
141 cout << "ir= " << ir << endl;
142 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
143 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
144 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
145 cout << "MAX====> " << maxprob << endl;
146 }
147 }
148 printf("MAXprob = %.10f\n",maxprob);
149*/
151 EvtVector4R D1 = p->getDaug(0)->getP4();
152 EvtVector4R D2 = p->getDaug(1)->getP4();
153 EvtVector4R D3 = p->getDaug(2)->getP4();
154
155 double P1[4], P2[4], P3[4];
156 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
157 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
158 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
159 //P1[0] = 0.840899; P1[1] = -0.624276; P1[2] = 0.0113601 ; P1[3] = -0.263898;
160 //P2[0] = 0.89303 ; P2[1] = 0.860432 ; P2[2] = -0.0491258; P2[3] = 0.187791 ;
161 //P3[0] = 0.288373; P3[1] = 0.102614 ; P3[2] = -0.132757 ; P3[3] = -0.191797;
162 //P1[0] = 0.817635; P1[1] = 0.140538 ; P1[2] = -0.501451 ; P1[3] = 0.386916;
163 //P2[0] = 0.952772; P2[1] = -0.164752; P2[2] = 0.437717 ; P2[3] = -0.818264;
164 //P3[0] = 0.260117; P3[1] = 0.213518 ; P3[2] = -0.0606403; P3[3] = -0.0132164;
165
166 double value;
167 int g0[5]={1,1,1,1,1};
168 int nstates=5;
169 calEva(P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value);
170
171 //value = calEva(P1, P2, P3);
172 setProb(value);
173 return ;
174}
175
176
177void EvtDsToKSpipi0::Com_Multi(double a1[2], double a2[2], double res[2])
178{
179 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
180 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
181}
182void EvtDsToKSpipi0::Com_Divide(double a1[2], double a2[2], double res[2])
183{
184 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
185 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
186 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
187}
188//------------base---------------------------------
189double EvtDsToKSpipi0::SCADot(double a1[4], double a2[4])
190{
191 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
192 return _cal;
193}
194double EvtDsToKSpipi0::barrier(int l, double sa, double sb, double sc, double r, double mass)
195{
196 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
197 if(q < 0) q = 1e-16;
198 double z;
199 z = q*r*r;
200 double sa0;
201 sa0 = mass*mass;
202 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
203 if(q0 < 0) q0 = 1e-16;
204 double z0 = q0*r*r;
205 double F = 0.0;
206 if(l == 0) F = 1;
207 if(l == 1) F = sqrt((1+z0)/(1+z));
208 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
209 return F;
210}
211
212void EvtDsToKSpipi0::calt1(double daug1[4], double daug2[4], double t1[4])
213{
214 double p, pq, tmp;
215 double pa[4], qa[4];
216 for(int i=0; i<4; i++) {
217 pa[i] = daug1[i] + daug2[i];
218 qa[i] = daug1[i] - daug2[i];
219 }
220 p = SCADot(pa,pa);
221 pq = SCADot(pa,qa);
222 tmp = pq/p;
223 for(int i=0; i<4; i++) {
224 t1[i] = qa[i] - tmp*pa[i];
225 }
226}
227void EvtDsToKSpipi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
228{
229 double p, r;
230 double pa[4], t1[4];
231 calt1(daug1,daug2,t1);
232 r = SCADot(t1,t1)/3.0;
233 for(int i=0; i<4; i++) {
234 pa[i] = daug1[i] + daug2[i];
235 }
236 p = SCADot(pa,pa);
237 for(int i=0; i<4; i++) {
238 for(int j=0; j<4; j++) {
239 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
240 }
241 }
242}
243//-------------------prop--------------------------------------------
244
245double EvtDsToKSpipi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
246{
247 double widm = 0.;
248 double m = sqrt(sa);
249 double tmp = sb-sc;
250 double tmp1 = sa+tmp;
251 double q = 0.25*tmp1*tmp1/sa-sb;
252 if(q<0) q = 1e-16;
253 double tmp2 = mass2+tmp;
254 double q0 = 0.25*tmp2*tmp2/mass2-sb;
255 if(q0<0) q0 = 1e-16;
256 double z = q*r2;
257 double z0 = q0*r2;
258 double t = q/q0;
259 if(l == 0) {widm = sqrt(t)*mass/m;}
260 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
261 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
262 return widm;
263}
264double EvtDsToKSpipi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
265{
266 double widm = 0.;
267 double m = sqrt(sa);
268 double tmp = sb-sc;
269 double tmp1 = sa+tmp;
270 double q = 0.25*tmp1*tmp1/sa-sb;
271 if(q<0) q = 1e-16;
272 double tmp2 = mass2+tmp;
273 double q0 = 0.25*tmp2*tmp2/mass2-sb;
274 if(q0<0) q0 = 1e-16;
275 double z = q*r2;
276 double z0 = q0*r2;
277 double F = (1+z0)/(1+z);
278 double t = q/q0;
279 widm = t*sqrt(t)*mass/m*F;
280 return widm;
281}
282void EvtDsToKSpipi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
283{
284 double a[2], b[2];
285 a[0] = 1;
286 a[1] = 0;
287 b[0] = mass2-sa;
288 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
289 Com_Divide(a,b,prop);
290}
291
292void EvtDsToKSpipi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
293{
294 double a[2], b[2];
295 double tmp = sb-sc;
296 double tmp1 = sa+tmp;
297 double q2 = 0.25*tmp1*tmp1/sa-sb;
298 if(q2<0) q2 = 1e-16;
299
300 double tmp2 = mass2+tmp;
301 double q02 = 0.25*tmp2*tmp2/mass2-sb;
302 if(q02<0) q02 = 1e-16;
303
304 double q = sqrt(q2);
305 double q0 = sqrt(q02);
306 double m = sqrt(sa);
307 double q03 = q0*q02;
308 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
309
310 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
311 double h0 = GS1*q0/mass*tmp3;
312 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
313 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
314 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
315
316 a[0] = 1.0+d*width/mass;
317 a[1] = 0.0;
318 b[0] = mass2-sa+width*f;
319 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
320 Com_Divide(a,b,prop);
321}
322
323
324void EvtDsToKSpipi0::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
325 const double m1430 = 1.441;
326 const double sa0 = 1.441*1.441; // m1430*m1430;
327 const double w1430 = 0.193;
328 const double Lass1 = 0.25/sa0;
329 double tmp = sb-sc;
330 double tmp1 = sa0+tmp;
331 double q0 = Lass1*tmp1*tmp1-sb;
332 if(q0<0) q0 = 1e-16;
333 double tmp2 = sa+tmp;
334 double qs = 0.25*tmp2*tmp2/sa-sb;
335 double q = sqrt(qs);
336 double width = w1430*q*m1430/sqrt(sa*q0);
337 double temp_R = atan(m1430*width/(sa0-sa));
338 if(temp_R<0) temp_R += math_pi;
339 double deltaR = -109.7 + temp_R;
340 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
341 if(temp_F<0) temp_F += math_pi;
342 double deltaF = 0.1 + temp_F;
343 double deltaS = deltaR + 2.0*deltaF;
344 double t1 = 0.96*sin(deltaF);
345 double t2 = sin(deltaR);
346 double CF[2], CS[2];
347 CF[0] = cos(deltaF);
348 CF[1] = sin(deltaF);
349 CS[0] = cos(deltaS);
350 CS[1] = sin(deltaS);
351 prop[0] = t1*CF[0] + t2*CS[0];
352 prop[1] = t1*CF[1] + t2*CS[1];
353}
354
355double EvtDsToKSpipi0::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
356 double pR[4], pD[4];
357 double temp_PDF, v_re;
358 temp_PDF = 0.0;
359 v_re = 0.0;
360 double B[2], s1, s2, s3, sR, sD;
361 for(int i=0; i<4; i++){
362 pR[i] = P1[i] + P2[i];
363 pD[i] = pR[i] + P3[i];
364 }
365 s1 = SCADot(P1,P1);
366 s2 = SCADot(P2,P2);
367 s3 = SCADot(P3,P3);
368 sR = SCADot(pR,pR);
369 sD = SCADot(pD,pD);
370 //int G[4][4];
371 //for(int i=0; i!=4; i++){
372 // for(int j=0; j!=4; j++){
373 // if(i==j){
374 // if(i==0) G[i][j] = 1;
375 // else G[i][j] = -1;
376 // }
377 // else G[i][j] = 0;
378 // }
379 //}
380 if(Ang == 0){
381 B[0] = 1;
382 B[1] = 1;
383 temp_PDF = 1;
384 }
385 if(Ang == 1){
386 B[0] = barrier(1,sR,s1,s2,3.0,mass);
387 B[1] = barrier(1,sD,sR,s3,5.0,1.9683);
388 double t1[4], T1[4];
389 calt1(P1,P2,t1);
390 calt1(pR,P3,T1);
391 temp_PDF = 0;
392 for(int i=0; i<4; i++){
393 temp_PDF += t1[i]*T1[i]*G[i][i];
394 }
395 }
396 if(Ang == 2){
397 B[0] = barrier(2,sR,s1,s2,3.0,mass);
398 B[1] = barrier(2,sD,sR,s3,5.0,1.9683);
399 double t2[4][4], T2[4][4];
400 calt2(P1,P2,t2);
401 calt2(pR,P3,T2);
402 temp_PDF = 0;
403 for(int i=0; i<4; i++){
404 for(int j=0; j<4; j++){
405 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
406 }
407 }
408 }
409 v_re = temp_PDF*B[0]*B[1];
410 return v_re;
411}
412
413
414void EvtDsToKSpipi0::calEva(double* Ks0, double* Pic, double* Pi0, double *mass1, double *width1, double *amp, double *phase,int* g0, int* modetype, int nstates, double & Result)
415{
416 //double Ks0[4], Pic[4], Pi0[4];
417 //double rRes = 9.0;
418 double P12[4], P23[4], P13[4];
419 double cof[2], amp_PDF[2], PDF[2];
420 double scpi, snpi, sks0;
421 double s12, s13, s23;
422 for(int i=0; i<4; i++){
423 P12[i] = Pic[i] + Ks0[i];
424 P13[i] = Pi0[i] + Ks0[i];
425 P23[i] = Pic[i] + Pi0[i];
426 }
427 scpi = SCADot(Pic,Pic);
428 snpi = SCADot(Pi0,Pi0);
429 sks0 = SCADot(Ks0,Ks0);
430 s12 = SCADot(P12,P12);
431 s13 = SCADot(P13,P13);
432 s23 = SCADot(P23,P23);
433 double pro[2], temp_PDF, amp_tmp[2];
434 double mass1sq;
435 amp_PDF[0] = 0;
436 amp_PDF[1] = 0;
437 PDF[0] = 0;
438 PDF[1] = 0;
439 amp_tmp[0] = 0;
440 amp_tmp[1] = 0;
441 for(int i=0; i<nstates; i++) {
442 amp_tmp[0] = 0;
443 amp_tmp[1] = 0;
444 mass1sq = mass1[i]*mass1[i];
445 cof[0] = amp[i]*cos(phase[i]);
446 cof[1] = amp[i]*sin(phase[i]);
447 temp_PDF = 0;
448 if(modetype[i] == 23||modetype[i] == 231){
449 temp_PDF = DDalitz( Pic, Pi0, Ks0, 1, mass1[i]);
450 if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],s23,scpi,snpi,9.0,pro);
451 if(g0[i]==0){
452 pro[0] = 1;
453 pro[1] = 0;
454 }
455 amp_tmp[0] = temp_PDF*pro[0];
456 amp_tmp[1] = temp_PDF*pro[1];
457
458 }
459 if(modetype[i] == 12){
460 temp_PDF = DDalitz(Ks0, Pic, Pi0, 1, mass1[i]);
461 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,sks0,scpi,rRes,1,pro);
462 if(g0[i]==12) KPiSLASS(s12,sks0,scpi,pro);
463 if(g0[i]==0){
464 pro[0] = 1;
465 pro[1] = 0;
466 }
467 amp_tmp[0] = temp_PDF*pro[0];
468 amp_tmp[1] = temp_PDF*pro[1];
469 }
470 if(modetype[i] == 13||modetype[i] == 131){
471 temp_PDF = DDalitz(Ks0, Pi0, Pic, 1, mass1[i]);
472 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,sks0,snpi,rRes,1,pro);
473 if(g0[i]==12) KPiSLASS(s13,sks0,snpi,pro);
474 if(g0[i]==0){
475 pro[0] = 1;
476 pro[1] = 0;
477 }
478 amp_tmp[0] = temp_PDF*pro[0];
479 amp_tmp[1] = temp_PDF*pro[1];
480 }
481 Com_Multi(amp_tmp,cof,amp_PDF);
482 PDF[0] += amp_PDF[0];
483 PDF[1] += amp_PDF[1];
484 }
485 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
486 if(value <=0) value = 1e-20;
487 Result = value;
488 //cout<<"value = "<<value<<" Result = "<<Result<<endl;
489}
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)
EvtDecayBase * clone()
virtual ~EvtDsToKSpipi0()
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