BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSKSpi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtDToKSKSpi.cc
11// the necessary file: EvtDToKSKSpi.cc
12//
13// Description: D+ -> KS0 KS0 pi+,
14//
15//
16// Modification history:
17// Liaoyuan Dong Thu Sep 29 22:06:58 CST 2022 Module created
18//------------------------------------------------------------------------
19//
20#include "EvtGenBase/EvtPatches.hh"
21#include "EvtGenBase/EvtParticle.hh"
22#include "EvtGenBase/EvtGenKine.hh"
23#include "EvtGenBase/EvtPDL.hh"
24#include "EvtGenBase/EvtReport.hh"
25#include "EvtGenModels/EvtDToKSKSpi.hh"
26#include "EvtGenBase/EvtComplex.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
28#include <stdlib.h>
29#include <iostream>
30#include <cmath>
31using namespace std;
32
33EvtDToKSKSpi::~EvtDToKSKSpi() {}
34
35void EvtDToKSKSpi::getName(std::string& model_name){
36 model_name="DToKSKSpi";
37}
38
42
43void EvtDToKSKSpi::init(){
44 checkNArg(0);
45 checkNDaug(3);
46 checkSpinParent(EvtSpinType::SCALAR);
47
48 mass[0] = 0.89167;//KpiSW
49 mass[1] = 0.89167;//K*892
50 width[0] = 0.05140;//KpiSW
51 width[1] = 0.05140;//K*892
52
53 rho[0] = 0.680648291453;//KpiSW
54 rho[1] = 1;//K*892
55 phi[0] = 1.837253259522;//KpiSW
56 phi[1] = 0;//K*892
57
58 spin[0] = 0; //KpiSW
59 spin[1] = 1; //K*892
60 modetype[0] = 13;//KpiSW
61 modetype[1] = 13;//K*892
62
63/*
64 std::cout << "EvtDToKSKSpi (May 01, 2023) ==> Initialization" << std::endl;
65 for (int i=0; i<2; i++) {
66 cout << i << "rho,phi = " << rho[i] << ", "<< phi[i] << endl;
67 }
68*/
69 math_pi = 3.1415926;
70 GS1 = 0.636619783;
71 GS2 = 0.01860182466;
72 GS3 = 0.1591549458; // 1/(2*math_2pi)
73 GS4 = 0.00620060822; // mass_Pion2/math_pi
74
75 int G[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
76
77
78}
79
80
81void EvtDToKSKSpi::initProbMax() {
82 setProbMax(361.0); // max 353.68 2022/11/08
83}
84
85
86
87void EvtDToKSKSpi::decay( EvtParticle *p){ //phsp
88
89//-------------------------------------------------begin-----------------------------------------
90/*
91//检查产生子正确后打开输出最大pdf
92 double maxprob = 0.0;
93 for(int ir=0;ir<=60000000;ir++){
94 p->initializePhaseSpace(getNDaug(),getDaugs());
95 EvtVector4R D1 = p->getDaug(0)->getP4();
96 EvtVector4R D2 = p->getDaug(1)->getP4();
97 EvtVector4R D3 = p->getDaug(2)->getP4();
98
99 double P1[4], P2[4], P3[4];
100 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
101 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
102 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
103
104 double value;
105 // value = calDalEva(P1, P2, P3);
106 int g0[2]={12,1};
107 int spin[2]={0,1};
108 int nstates=2;
109 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
110 if (value<0) continue;
111 if(value>maxprob) {
112 maxprob=value;
113 cout << "ir= " << ir << endl;
114 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
115 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
116 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
117 cout << "MAX====> " << maxprob << endl;
118 }
119 }
120 printf("MAXprob = %.10f\n",maxprob);
121*/
122//---------------------------------------------endls-----------------------------------------------------
123
124 //检验最大值关闭
125 p->initializePhaseSpace(getNDaug(),getDaugs());
126 EvtVector4R D1 = p->getDaug(0)->getP4();
127 EvtVector4R D2 = p->getDaug(1)->getP4();
128 EvtVector4R D3 = p->getDaug(2)->getP4();
129
130 double P1[4], P2[4], P3[4];
131 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
132 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
133 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
134 //Check
135 //P1[0] = 0.817527;P1[1] = -0.0302783; P1[2] = 0.535352; P1[3] = -0.364866;
136 //P2[0] = 0.660657;P2[1] = 0.132324; P2[2] = -0.368085; P2[3] = 0.18912 ;
137 //P3[0] = 0.40678; P3[1] = 0.0991563 ; P3[2] = -0.281318; P3[3] = 0.238786;
138
139 double value;
140 int g0[2]={12,1};
141 int spin[2]={0,1};
142 int nstates=2;
143 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
144 setProb(value); //检验最大值关闭
145
146 return ;
147}
148
149void EvtDToKSKSpi::Com_Multi(double a1[2], double a2[2], double res[2])
150{
151 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
152 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
153}
154void EvtDToKSKSpi::Com_Divide(double a1[2], double a2[2], double res[2])
155{
156 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
157 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
158 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
159}
160//------------base---------------------------------
161double EvtDToKSKSpi::SCADot(double a1[4], double a2[4])
162{
163 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
164 return _cal;
165}
166double EvtDToKSKSpi::barrier(int l, double sa, double sb, double sc, double r, double mass)
167{
168 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
169 if(q < 0) q = 1e-16;
170 double z;
171 z = q*r*r;
172 double sa0;
173 sa0 = mass*mass;
174 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
175 if(q0 < 0) q0 = 1e-16;
176 double z0 = q0*r*r;
177 double F = 0.0;
178 if(l == 0) F = 1;
179 if(l == 1) F = sqrt((1+z0)/(1+z));
180 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
181 return F;
182}
183
184void EvtDToKSKSpi::calt1(double daug1[4], double daug2[4], double t1[4])
185{
186 double p, pq, tmp;
187 double pa[4], qa[4];
188 for(int i=0; i<4; i++) {
189 pa[i] = daug1[i] + daug2[i];
190 qa[i] = daug1[i] - daug2[i];
191 }
192 p = SCADot(pa,pa);
193 pq = SCADot(pa,qa);
194 tmp = pq/p;
195 for(int i=0; i<4; i++) {
196 t1[i] = qa[i] - tmp*pa[i];
197 }
198}
199
200void EvtDToKSKSpi::calt2(double daug1[4], double daug2[4], double t2[4][4])
201{
202 double p, r;
203 double pa[4], t1[4];
204 calt1(daug1,daug2,t1);
205 r = SCADot(t1,t1)/3.0;
206 for(int i=0; i<4; i++) {
207 pa[i] = daug1[i] + daug2[i];
208 }
209 p = SCADot(pa,pa);
210 for(int i=0; i<4; i++) {
211 for(int j=0; j<4; j++) {
212 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
213 }
214 }
215}
216
217double EvtDToKSKSpi::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
218{
219 double widm = 0.;
220 double m = sqrt(sa);
221 double tmp = sb-sc;
222 double tmp1 = sa+tmp;
223 double q = 0.25*tmp1*tmp1/sa-sb;
224 if(q<0) q = 1e-16;
225 double tmp2 = mass2+tmp;
226 double q0 = 0.25*tmp2*tmp2/mass2-sb;
227 if(q0<0) q0 = 1e-16;
228 double z = q*r2;
229 double z0 = q0*r2;
230 double t = q/q0;
231 if(l == 0) {widm = sqrt(t)*mass/m;}
232 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
233 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
234 return widm;
235}
236double EvtDToKSKSpi::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
237{
238 double widm = 0.;
239 double m = sqrt(sa);
240 double tmp = sb-sc;
241 double tmp1 = sa+tmp;
242 double q = 0.25*tmp1*tmp1/sa-sb;
243 if(q<0) q = 1e-16;
244 double tmp2 = mass2+tmp;
245 double q0 = 0.25*tmp2*tmp2/mass2-sb;
246 if(q0<0) q0 = 1e-16;
247 double z = q*r2;
248 double z0 = q0*r2;
249 double F = (1+z0)/(1+z);
250 double t = q/q0;
251 widm = t*sqrt(t)*mass/m*F;
252 return widm;
253}
254void EvtDToKSKSpi::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
255{
256 double a[2], b[2];
257 a[0] = 1;
258 a[1] = 0;
259 b[0] = mass2-sa;
260 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
261 Com_Divide(a,b,prop);
262}
263
264void EvtDToKSKSpi::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
265{
266 double a[2], b[2];
267 double tmp = sb-sc;
268 double tmp1 = sa+tmp;
269 double q2 = 0.25*tmp1*tmp1/sa-sb;
270 if(q2<0) q2 = 1e-16;
271
272 double tmp2 = mass2+tmp;
273 double q02 = 0.25*tmp2*tmp2/mass2-sb;
274 if(q02<0) q02 = 1e-16;
275
276 double q = sqrt(q2);
277 double q0 = sqrt(q02);
278 double m = sqrt(sa);
279 double q03 = q0*q02;
280 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
281
282 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
283 double h0 = GS1*q0/mass*tmp3;
284 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
285 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
286 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
287
288 a[0] = 1.0+d*width/mass;
289 a[1] = 0.0;
290 b[0] = mass2-sa+width*f;
291 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
292 Com_Divide(a,b,prop);
293}
294
295
296void EvtDToKSKSpi::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
297 const double m1430 = 1.441;
298 const double sa0 = 1.441*1.441; // m1430*m1430;
299 const double w1430 = 0.193;
300 const double Lass1 = 0.25/sa0;
301 double tmp = sb-sc;
302 double tmp1 = sa0+tmp;
303 double q0 = Lass1*tmp1*tmp1-sb;
304 //if(q0<0) q0 = 1e-16;
305 double tmp2 = sa+tmp;
306 double qs = 0.25*tmp2*tmp2/sa-sb;
307 double q = sqrt(qs);
308 double width = w1430*q*m1430/sqrt(sa*q0);
309 double temp_R = atan(m1430*width/(sa0-sa));
310 if(temp_R<0) temp_R += math_pi;
311 double deltaR = -109.7 * math_pi / 180 + temp_R;
312 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
313 if(temp_F<0) temp_F += math_pi;
314 double deltaF = 0.1 * math_pi / 180 + temp_F;
315 double deltaS = deltaR + 2.0*deltaF;
316 double t1 = 0.96*sin(deltaF);
317 double t2 = sin(deltaR);
318 double CF[2], CS[2];
319 CF[0] = cos(deltaF);
320 CF[1] = sin(deltaF);
321 CS[0] = cos(deltaS);
322 CS[1] = sin(deltaS);
323 prop[0] = t1*CF[0] + t2*CS[0];
324 prop[1] = t1*CF[1] + t2*CS[1];
325}
326
327double EvtDToKSKSpi::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
328 double pR[4], pD[4];
329 double temp_PDF, v_re;
330 temp_PDF = 0.0;
331 v_re = 0.0;
332 double B[2], s1, s2, s3, sR, sD;
333 for(int i=0; i<4; i++){
334 pR[i] = P1[i] + P2[i];
335 pD[i] = pR[i] + P3[i];
336 }
337 s1 = SCADot(P1,P1);
338 s2 = SCADot(P2,P2);
339 s3 = SCADot(P3,P3);
340 sR = SCADot(pR,pR);
341 sD = SCADot(pD,pD);
342
343 int G[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
344
345 if(Ang == 0){
346 B[0] = 1;
347 B[1] = 1;
348 temp_PDF = 1;
349 }
350 if(Ang == 1){
351 B[0] = barrier(1,sR,s1,s2,3.0,mass);
352 B[1] = barrier(1,sD,sR,s3,5.0,1.86966);
353 double t1[4], T1[4];
354 calt1(P1,P2,t1);
355 calt1(pR,P3,T1);
356 temp_PDF = 0;
357 for(int i=0; i<4; i++){
358 temp_PDF += t1[i]*T1[i]*G[i][i];
359 }
360 }
361 if(Ang == 2){
362 B[0] = barrier(2,sR,s1,s2,3.0,mass);
363 B[1] = barrier(2,sD,sR,s3,5.0,1.86966);
364 double t2[4][4], T2[4][4];
365 calt2(P1,P2,t2);
366 calt2(pR,P3,T2);
367 temp_PDF = 0;
368 for(int i=0; i<4; i++){
369 for(int j=0; j<4; j++){
370 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
371 }
372 }
373 }
374 v_re = temp_PDF*B[0]*B[1];
375 return v_re;
376}
377
378void EvtDToKSKSpi::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)
379{
380 double rRes = 9.0;
381 double P12[4], P23[4], P13[4];
382 double cof[2], amp_PDF[2], PDF[2];
383 double scpi, sks02, sks01;
384 double s12, s13, s23;
385 for(int i=0; i<4; i++){
386 P12[i] = Ks01[i] + Ks02[i];
387 P13[i] = Ks01[i] + Pic[i];
388 P23[i] = Ks02[i] + Pic[i];
389 }
390 sks01 = SCADot(Ks01,Ks01);
391 sks02 = SCADot(Ks02,Ks02);
392 scpi = SCADot(Pic,Pic);
393 s12 = SCADot(P12,P12);
394 s13 = SCADot(P13,P13);
395 s23 = SCADot(P23,P23);
396 double pro[2], temp_PDF, amp_tmp[2],temp_PDF1 ,temp_PDF2,pro1[2],pro2[2];
397 double mass1sq;
398 amp_PDF[0] = 0;
399 amp_PDF[1] = 0;
400 PDF[0] = 0;
401 PDF[1] = 0;
402 amp_tmp[0] = 0;
403 amp_tmp[1] = 0;
404 for(int i=0; i<nstates; i++) {
405 amp_tmp[0] = 0;
406 amp_tmp[1] = 0;
407 mass1sq = mass1[i]*mass1[i];
408 cof[0] = amp[i]*cos(phase[i]);
409 cof[1] = amp[i]*sin(phase[i]);
410 temp_PDF = 0;
411 if(modetype[i] == 12){
412 temp_PDF = DDalitz(Ks01, Ks02, Pic, spin[i], mass1[i]);
413 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,sks01,sks02,rRes,spin[i],pro);
414 if(g0[i]==12) KPiSLASS(s12,sks01,sks02,pro);
415 if(g0[i]==0){
416 pro[0] = 1;
417 pro[1] = 0;
418 }
419 amp_tmp[0] = temp_PDF*pro[0];
420 amp_tmp[1] = temp_PDF*pro[1];
421 }
422 if(modetype[i] == 13){
423 temp_PDF1 = DDalitz(Ks01, Pic, Ks02, spin[i], mass1[i]);
424 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,sks01,scpi,rRes,spin[i],pro1);
425 if(g0[i]==12) KPiSLASS(s13,sks01,scpi,pro1);
426 if(g0[i]==0){
427 pro1[0] = 1;
428 pro1[1] = 0;
429 }
430 temp_PDF2 = DDalitz(Ks02, Pic, Ks01, spin[i], mass1[i]);
431 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s23,sks02,scpi,rRes,spin[i],pro2);
432 if(g0[i]==12) KPiSLASS(s23,sks02,scpi,pro2);
433 if(g0[i]==0){
434 pro2[0] = 1;
435 pro2[1] = 0;
436 }
437 amp_tmp[0] = ( temp_PDF1*pro1[0]+temp_PDF2*pro2[0]);
438 amp_tmp[1] = ( temp_PDF1*pro1[1]+temp_PDF2*pro2[1]);
439 }
440 Com_Multi(amp_tmp,cof,amp_PDF);
441 PDF[0] += amp_PDF[0];
442 PDF[1] += amp_PDF[1];
443
444 }
445 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
446 if(value <=0) value = 1e-20;
447 Result = value;
448//Check
449 //printf("%s %8.15f\n","AAAvalue = ",value);
450}
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
EvtDToKSKSpi::getName
void getName(std::string &name)
Definition EvtDToKSKSpi.cc:35
EvtDToKSKSpi::decay
void decay(EvtParticle *p)
Definition EvtDToKSKSpi.cc:87
EvtDToKSKSpi::~EvtDToKSKSpi
virtual ~EvtDToKSKSpi()
Definition EvtDToKSKSpi.cc:33
EvtDToKSKSpi::clone
EvtDecayBase * clone()
Definition EvtDToKSKSpi.cc:39
EvtDToKSKSpi::initProbMax
void initProbMax()
Definition EvtDToKSKSpi.cc:81
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