BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKpiEtap.cc
Go to the documentation of this file.
1// Environment:
2// This software is part of the EvtGen package developed jointly
3// for the BaBar and CLEO collaborations. If you use all or part
4// of it, please give an appropriate acknowledgement.
5//
6// Copyright Information: See EvtGen/COPYRIGHT
7// Copyright (C) 1998 Caltech, UCSB
8//
9// Module: EvtD0ToEtapipi.cc
10//
11// Description: Routine to handle three-body decays of D0/D0_bar or D+/D-
12//
13// Modification history:
14//
15// Liaoyuan Dong Jul 3 23:56:42 2023 Module created
16//
17//------------------------------------------------------------------------
18//
19#include "EvtGenBase/EvtPatches.hh"
20#include "EvtGenBase/EvtParticle.hh"
21#include "EvtGenBase/EvtGenKine.hh"
22#include "EvtGenBase/EvtPDL.hh"
23#include "EvtGenBase/EvtReport.hh"
24#include "EvtGenBase/EvtVector4R.hh"
25#include "EvtGenBase/EvtComplex.hh"
26#include "EvtGenModels/EvtD0ToKpiEtap.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
28#include <stdlib.h>
29using namespace std;
30
31EvtD0ToKpiEtap::~EvtD0ToKpiEtap() {}
32
33void EvtD0ToKpiEtap::getName(std::string& model_name){
34 model_name="D0ToKpiEtap";
35}
36
40void EvtD0ToKpiEtap::init(){
41 checkNArg(0);
42 checkNDaug(3);
43 checkSpinParent(EvtSpinType::SCALAR);
44 checkSpinDaughter(0,EvtSpinType::SCALAR);
45 checkSpinDaughter(1,EvtSpinType::SCALAR);
46 checkSpinDaughter(2,EvtSpinType::SCALAR);
47
48 phi[0] = 0.0;
49 rho[0] = 1.0;
50 phi[1] = 3.7295;
51 rho[1] = 2.2532e-01; //892
52 phi[2] = 0.6890853;
53 rho[2] = 1.6780e+00;//1680
54 modetype[0]= 132;
55 modetype[1]= 12;
56 modetype[2]= 12;
57
58 mass[2] = 1.718;
59 mass[1] = 0.89555;
60 width[2]= 0.322;
61 width[1]= 0.0473;
62 mass[0] = 1.414;
63 width[0]= 0.232;
64 const double mk0 = 0.497614;
65 const double mass_Kaon = 0.49368;
66 const double mass_Pion = 0.13957;
67 const double mass_Pi0 = 0.1349766;
68 const double meta = 0.547862;
69 const double metap = 0.95778;
70 mpi = 0.13957;
71 mD = 1.86486;
72 mDs = 1.9683;
73 rRes= 9.0;
74 rD = 5.0;
75 math_pi = 3.1415926;
76 GS1 = 0.636619783;
77 GS2 = 0.01860182466;
78 GS3 = 0.1591549458; // 1/(2*math_2pi)
79 GS4 = 0.00620060822; // mass_Pion2/math_pi
80
81 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
82 for (int i=0; i<4; i++) {
83 for (int j=0; j<4; j++) {
84 G[i][j] = GG[i][j];
85 }
86 }
87}
91
92void EvtD0ToKpiEtap::decay( EvtParticle *p){
93
94 // This piece of code could in principle be used to calculate maximum
95 // probablity on fly. But as it uses high number of points and model
96 // deals with single final state, we keep hardcoded number for now rather
97 // than adapting code to work here.
98/*
99 double maxprob = 0.0;
100 for(int ir=0;ir<=60000000;ir++){
101 p->initializePhaseSpace(getNDaug(),getDaugs());
102 EvtVector4R D1 = p->getDaug(0)->getP4();
103 EvtVector4R D2 = p->getDaug(1)->getP4();
104 EvtVector4R D3 = p->getDaug(2)->getP4();
105
106 double P1[4], P2[4], P3[4];
107 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
108 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
109 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
110
111 double value;
112 int g0[3]={0,1,1};
113 int spin[3]={0,1,1};
114 int nstates=3;
115 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
116 if(value>maxprob) {
117 maxprob=value;
118 cout << "ir = " << ir << " maxProb= " << value << endl;
119 }
120 }
121 cout << "maxProb = " << maxprob << endl;
122*/
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 int g0[3]={0,1,1};
135 int spin[3]={0,1,1};
136 int nstates=3;
137 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
138
139 setProb(value);
140 return ;
141}
142
143void EvtD0ToKpiEtap::Com_Multi(double a1[2], double a2[2], double res[2])
144{
145 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
146 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
147}
148void EvtD0ToKpiEtap::Com_Divide(double a1[2], double a2[2], double res[2])
149{
150 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
151 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
152 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
153}
154//------------base---------------------------------
155double EvtD0ToKpiEtap::SCADot(double a1[4], double a2[4])
156{
157 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
158 return _cal;
159}
160double EvtD0ToKpiEtap::barrier(int l, double sa, double sb, double sc, double r, double mass)
161{
162 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
163 if(q < 0) q = -q;
164 double z;
165 z = q*r*r;
166 double sa0;
167 sa0 = mass*mass;
168 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
169 if(q0 < 0) q0 = -1*q0;
170 double z0 = q0*r*r;
171 double F = 0.0;
172 if(l == 0) F = 1;
173 if(l == 1) F = sqrt((1+z0)/(1+z));
174 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
175 return F;
176}
177void EvtD0ToKpiEtap::calt1(double daug1[4], double daug2[4], double t1[4])
178{
179 double p, pq, tmp;
180 double pa[4], qa[4];
181 for(int i=0; i<4; i++) {
182 pa[i] = daug1[i] + daug2[i];
183 qa[i] = daug1[i] - daug2[i];
184 }
185 p = SCADot(pa,pa);
186 pq = SCADot(pa,qa);
187 tmp = pq/p;
188 for(int i=0; i<4; i++) {
189 t1[i] = qa[i] - tmp*pa[i];
190 }
191}
192void EvtD0ToKpiEtap::calt2(double daug1[4], double daug2[4], double t2[4][4])
193{
194 double p, r;
195 double pa[4], t1[4];
196 calt1(daug1,daug2,t1);
197 r = SCADot(t1,t1)/3.0;
198 for(int i=0; i<4; i++) {
199 pa[i] = daug1[i] + daug2[i];
200 }
201 p = SCADot(pa,pa);
202 for(int i=0; i<4; i++) {
203 for(int j=0; j<4; j++) {
204 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
205 }
206 }
207}
208//-------------------prop--------------------------------------------
209double EvtD0ToKpiEtap::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
210{
211 double widm = 0.;
212 double m = sqrt(sa);
213 double tmp = sb-sc;
214 double tmp1 = sa+tmp;
215 double q = 0.25*tmp1*tmp1/sa-sb;
216 if(q<0) q = -q;
217 double tmp2 = mass2+tmp;
218 double q0 = 0.25*tmp2*tmp2/mass2-sb;
219 if(q0<0) q0 = -q0;
220 double z = q*r2;
221 double z0 = q0*r2;
222 double t = q/q0;
223 if(l == 0) {widm = sqrt(t)*mass/m;}
224 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
225 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
226 return widm;
227}
228double EvtD0ToKpiEtap::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
229{
230 double widm = 0.;
231 double m = sqrt(sa);
232 double tmp = sb-sc;
233 double tmp1 = sa+tmp;
234 double q = 0.25*tmp1*tmp1/sa-sb;
235 if(q<0) q = -q;
236 double tmp2 = mass2+tmp;
237 double q0 = 0.25*tmp2*tmp2/mass2-sb;
238 if(q0<0) q0 = -q0;
239 double z = q*r2;
240 double z0 = q0*r2;
241 double F = (1+z0)/(1+z);
242 double t = q/q0;
243 widm = t*sqrt(t)*mass/m*F;
244 return widm;
245}
246void EvtD0ToKpiEtap::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
247{
248 double a[2], b[2];
249 a[0] = 1;
250 a[1] = 0;
251 double q=0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
252
253 b[0] = mass2-sa;
254 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
255 Com_Divide(a,b,prop);
256}
257
258void EvtD0ToKpiEtap::propagatorFlatte(double mass, double width, double sa, double sb, double sc, int r, double prop[2]){
259 double q, qKsK,qetapi;
260 // double qKsK,qetapi;
261 double rhoab[2], rhoKsK[2];
262 q = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
263 qetapi=0.25*(sa+0.547862-0.13957039)*(sa+0.547862-0.13957039)/sa-0.547862*0.547862;
264 if(r == 0) qKsK = 0.25*sa - 0.49368*0.49368;
265 if(r == 1) qKsK = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa - sb;
266 if(qetapi>0){
267 rhoab[0] = 2*sqrt(qetapi/sa);
268 rhoab[1] = 0;
269 }
270 if(qetapi<0){
271 rhoab[0] = 0;
272 rhoab[1] = 2*sqrt(-qetapi/sa);
273 }
274 if(qKsK>0){
275 rhoKsK[0] = 2*sqrt(qKsK/sa);
276 rhoKsK[1] = 0;
277 }
278 if(qKsK<0){
279 rhoKsK[0] = 0;
280 rhoKsK[1] = 2*sqrt(-qKsK/sa);
281 }
282 double a[2], b[2];
283 a[0] = 1;
284 a[1] = 0;
285 b[0] = mass*mass - sa + 0.341*rhoab[1] + 0.892*0.341*rhoKsK[1];
286 b[1] = - (0.341*rhoab[0] + 0.892*0.341*rhoKsK[0]);
287 Com_Divide(a,b,prop);
288}
289
290
291void EvtD0ToKpiEtap::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
292{
293 double a[2], b[2];
294 double tmp = sb-sc;
295 double tmp1 = sa+tmp;
296 double q2 = 0.25*tmp1*tmp1/sa-sb;
297 if(q2<0) q2 = 1e-16;
298
299 double tmp2 = mass2+tmp;
300 double q02 = 0.25*tmp2*tmp2/mass2-sb;
301 if(q02<0) q02 = 1e-16;
302
303 double q = sqrt(q2);
304 double q0 = sqrt(q02);
305 double m = sqrt(sa);
306 double q03 = q0*q02;
307 //double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
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
323void EvtD0ToKpiEtap::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
324 const double m1430 = 1.441;
325 const double sa0 = 1.441*1.441; // m1430*m1430;
326 const double w1430 = 0.193;
327 const double Lass1 = 0.25/sa0;
328 double tmp = sb-sc;
329 double tmp1 = sa0+tmp;
330 double q0 = Lass1*tmp1*tmp1-sb;
331 if(q0<0) q0 = -1*q0;
332 double tmp2 = sa+tmp;
333 double qs = 0.25*tmp2*tmp2/sa-sb;
334 double q = sqrt(qs);
335 double width = w1430*q*m1430/sqrt(sa*q0);
336 double temp_R = atan(m1430*width/(sa0-sa));
337 if(temp_R<0) temp_R += math_pi;
338 double deltaR = -1.915 + temp_R; //fiR=-109.7
339 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
340 if(temp_F<0) temp_F += math_pi;
341 double deltaF = 0.002 + temp_F; //fiF=0.1
342 double deltaS = deltaR + 2.0*deltaF;
343 double t1 = 0.96*sin(deltaF);
344 double t2 = sin(deltaR);
345 double CF[2], CS[2];
346 CF[0] = cos(deltaF);
347 CF[1] = sin(deltaF);
348 CS[0] = cos(deltaS);
349 CS[1] = sin(deltaS);
350 prop[0] = t1*CF[0] + t2*CS[0];
351 prop[1] = t1*CF[1] + t2*CS[1];
352}
353double EvtD0ToKpiEtap::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
354 double pR[4], pD[4];
355 double temp_PDF, v_re;
356 temp_PDF = 0.0;
357 v_re = 0.0;
358 double B[2], s1, s2, s3, sR, sD;
359 for(int i=0; i<4; i++){
360 pR[i] = P1[i] + P2[i];
361 pD[i] = pR[i] + P3[i];
362 }
363 s1 = SCADot(P1,P1);
364 s2 = SCADot(P2,P2);
365 s3 = SCADot(P3,P3);
366 sR = SCADot(pR,pR);
367 sD = SCADot(pD,pD);
368 if(Ang == 0){
369 B[0] = 1;
370 B[1] = 1;
371 temp_PDF = 1;
372 }
373 if(Ang == 1){
374 B[0] = barrier(1,sR,s1,s2,3.0,mass);
375 B[1] = barrier(1,sD,sR,s3,5.0,1.86484);
376 //B[0] = Barrier(1,sR,s1,s2,9.0);
377 //B[1] = Barrier(1,sD,sR,s3,25.0);
378 double t1[4], T1[4];
379 calt1(P1,P2,t1);
380 calt1(pR,P3,T1);
381 temp_PDF = 0;
382 for(int i=0; i<4; i++){
383 temp_PDF += t1[i]*T1[i]*G[i][i];
384 }
385 }
386 if(Ang == 2){
387 B[0] = barrier(2,sR,s1,s2,3.0,mass);
388 B[1] = barrier(2,sD,sR,s3,5.0,1.86484);
389 // B[0] = Barrier(2,sR,s1,s2,9.0);
390 // B[1] = Barrier(2,sD,sR,s3,25.0);
391 double t2[4][4], T2[4][4];
392 calt2(P1,P2,t2);
393 calt2(pR,P3,T2);
394 temp_PDF = 0;
395 for(int i=0; i<4; i++){
396 for(int j=0; j<4; j++){
397 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
398 }
399 }
400 }
401 v_re = temp_PDF*B[0]*B[1];
402 return v_re;
403}
404
405
406
407void EvtD0ToKpiEtap::calEva(double* Ks0, double* Kc, double* Pi0, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result)
408{
409 double P12[4], P23[4], P13[4];
410 double cof[2], amp_PDF[2], PDF[2];
411 double snpi, sck, sks0;
412 double s12, s13, s23;
413 for(int i=0; i<4; i++){
414 P12[i] = Kc[i] + Ks0[i];
415 P13[i] = Pi0[i] + Ks0[i];
416 P23[i] = Kc[i] + Pi0[i];
417 }
418 sck = SCADot(Kc,Kc);
419 snpi = SCADot(Pi0,Pi0);
420 sks0 = SCADot(Ks0,Ks0);
421 s12 = SCADot(P12,P12);
422 s13 = SCADot(P13,P13);
423 s23 = SCADot(P23,P23);
424 double pro[2], temp_PDF, amp_tmp[2],temp_PDF1 ,temp_PDF2,pro1[2],pro2[2];
425 double mass1sq;
426 amp_PDF[0] = 0;
427 amp_PDF[1] = 0;
428 PDF[0] = 0;
429 PDF[1] = 0;
430 amp_tmp[0] = 0;
431 amp_tmp[1] = 0;
432 for(int i=0; i<nstates; i++) {
433 amp_tmp[0] = 0;
434 amp_tmp[1] = 0;
435 mass1sq = mass1[i]*mass1[i];
436 cof[0] = amp[i]*cos(phase[i]);
437 cof[1] = amp[i]*sin(phase[i]);
438 temp_PDF = 0;
439
440 if(modetype[i] == 12){
441 temp_PDF = DDalitz(Ks0, Kc, Pi0, spin[i], mass1[i]);
442 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,sks0,sck,rRes,spin[i],pro);
443 // if(g0[i]==2) KPiSLASS(s12,sks0,scpi,pro);
444 if(g0[i]==2) {propagatorFlatte(mass1[i],width1[i],s12,sks0,sck,1,pro);
445 pro[0]=pro[0]*(0.01+0.990*0.990)/(0.01+s12);
446 pro[1]=pro[1]*(0.01+0.990*0.990)/(0.01+s12);
447 // pro[0]=pro[0]*(0.6*0.6)/(0.6*0.6)+(s12-0.990*0.990)*(s12-0.990*0.990);
448 // pro[1]=pro[1]*(0.6*0.6)/(0.6*0.6)+(s12-0.990*0.990)*(s12-0.990*0.990);
449 }
450 if(g0[i]==3) propagatorFlatte(mass1[i],width1[i],s12,sks0,sck,1,pro);//Only for a0(980)
451
452 if(g0[i]==0){
453 pro[0] = 1;
454 pro[1] = 0;
455 }
456 amp_tmp[0] = temp_PDF*pro[0];
457 amp_tmp[1] = temp_PDF*pro[1];
458 }
459
460 if(modetype[i] == 132){
461 KPiSLASS(s12,sks0,sck,pro);
462 amp_tmp[0] = pro[0];
463 amp_tmp[1] = pro[1];
464
465
466
467 }
468 Com_Multi(amp_tmp,cof,amp_PDF);
469 PDF[0] += amp_PDF[0];
470 PDF[1] += amp_PDF[1];
471
472 }
473 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
474 if(value <=0) value = 1e-20;
475 Result = value;
476}
477
478
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")
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
EvtD0ToKpiEtap::decay
void decay(EvtParticle *p)
Definition EvtD0ToKpiEtap.cc:92
EvtD0ToKpiEtap::clone
EvtDecayBase * clone()
Definition EvtD0ToKpiEtap.cc:37
EvtD0ToKpiEtap::~EvtD0ToKpiEtap
virtual ~EvtD0ToKpiEtap()
Definition EvtD0ToKpiEtap.cc:31
EvtD0ToKpiEtap::getName
void getName(std::string &name)
Definition EvtD0ToKpiEtap.cc:33
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:533
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:520
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:297
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:504
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:65
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:482
EvtDecayProb::setProb
void setProb(double prob)
Definition EvtDecayProb.hh:34
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:120
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:84
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1070
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:179
s12
double double double double * s12
Definition qcdloop1.h:77
s23
double double double double double * s23
Definition qcdloop1.h:77
b
const double b
Definition slope.cxx:9