BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKSpi0pi0.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: EvtD0ToKSpi0pi0.cc
11//
12// Description: D0 -> KS pi0 pi0
13//
14// Modification history:
15//
16// Liaoyuan Dong Feb 19 22:49:36 2023 Module created
17// Aug 22 21:58:51 2023 Use pipi S-wave K-matrix
18//
19//------------------------------------------------------------------------
21#include <fstream>
22#include <stdlib.h>
25#include "EvtGenBase/EvtPDL.hh"
30#include <string>
31#include <vector>
32#include <math.h>
33#include "TMath.h"
36using std::endl;
37
39
40void EvtD0ToKSpi0pi0::getName(std::string& model_name){
41 model_name="D0ToKSpi0pi0";
42}
43
47
49 // check that there are 0 arguments
50 checkNArg(0);
51 checkNDaug(3);
56
57 pi180inv = 1.0/EvtConst::radToDegrees;
58 phi[0] = 0.0;
59 phi[1] = -1.1735907962948264327;
60 phi[2] = -0.98738790124871123055;
61 phi[3] = 2.4055643732573601667;
62 phi[4] = 1.9872186812020204982;
63 phi[5] = 0.25034666009967576628;
64 phi[6] = -2.4035171276352995662;
65
66 rho[0] = 1.0;
67 rho[1] = 1.407574144549414541;
68 rho[2] = 1.3964570662861071071;
69 rho[3] = 2.5234814225240658203;
70 rho[4] = 2.0222401273514734044;
71 rho[5] = 0.69741141764274061643;
72 rho[6] = 0.84723998937758082661;
73
74 modetype[0]= 12;
75 modetype[1]= 12;
76 modetype[2]= 12;
77 modetype[3]= 12;
78 modetype[4]= 12;
79 modetype[5]= 23;
80 modetype[6]= 23;
81
82 //std::cout << "Initializing EvtD0ToKSpi0pi0" << std::endl;
83 //for (int i=0; i<7; i++) {
84 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
85 //}
86
87 width[0] = 0.047300000000000001765;
88 width[1] = 0.23200000000000001177;
89 width[2] = 0.10899999999999999967;
90 width[3] = 0.0084899999999999992834;
91 width[4] = 0.32200000000000000844;
92 width[5] = 0.020000000000000000416;
93 width[6] = 0.18670000000000000484;
94
95 mass[0] = 0.89554999999999995719;
96 mass[1] = 1.4139999999999999236;
97 mass[2] = 1.4323999999999998956;
98 mass[3] = 1.3999999999999999112;
99 mass[4] = 1.7179999999999999716;
100 mass[5] = 0.9000000000000000222;
101 mass[6] = 1.2755000000000000782;
102
103 mD0 = 1.86483;
104 mK0 = 0.497611;
105 mKa = 0.49368;
106 mPi = 0.13957;
107 meta = 0.54775;
108 mK02 = 0.237616707;
109 mPi2 = 0.01947978;
110 mass_EtaP = 0.95778;
111 mass_Kaon = 0.49368;
112
113 math_pi = 3.1415926;
114 mass_Pion = 0.13957;
115 mass_Pion2 = 0.0194797849;
116 mass_2Pion = 0.27914;
117 math_2pi = 6.2831852;
118 rD2 = 25.0; // 5*5
119 rRes2 = 9.0; // 3*3
120 g1 = 0.5468;
121 g2 = 0.23;
122 GS1 = 0.636619783;
123 GS2 = 0.01860182466;
124 GS3 = 0.1591549458;
125 GS4 = 0.00620060822;
126
127 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
128 for (int i=0; i<4; i++) {
129 for (int j=0; j<4; j++) {
130 G[i][j] = GG[i][j];
131 }
132 }
133}
134
136 setProbMax(944.0);//MAXprob = 943.6824670811
137}
138
140/*
141 double maxprob = 0.0;
142 for(int ir=0;ir<=60000000;ir++){
143 p->initializePhaseSpace(getNDaug(),getDaugs());
144 EvtVector4R D1 = p->getDaug(0)->getP4();
145 EvtVector4R D2 = p->getDaug(1)->getP4();
146 EvtVector4R D3 = p->getDaug(2)->getP4();
147
148 double P1[4], P2[4], P3[4];
149 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
150 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
151 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
152
153 double value;
154 int nstates=7;
155 int g0[7]={1,1,1,4,1,6,2};
156 int spin[7]={1,1,2,0,1,0,2};
157 double r0[7]={3.0,3.0,3.0,3.0,3.0,3.0,3.0};
158 double r1[7]={5.0,5.0,5.0,5.0,5.0,5.0,5.0};
159
160 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
161
162 if (value<0) continue;
163 if(value>maxprob) {
164 maxprob=value;
165 cout << "ir= " << ir << endl;
166 cout << "MAX====> " << maxprob << endl;
167 }
168 }
169 printf("MAXprob = %.10f\n",maxprob);
170*/
172 EvtVector4R D1 = p->getDaug(0)->getP4();
173 EvtVector4R D2 = p->getDaug(1)->getP4();
174 EvtVector4R D3 = p->getDaug(2)->getP4();
175
176 double P1[4], P2[4], P3[4];
177 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
178 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
179 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
180
181 double value;
182 int nstates=7;
183 int g0[7]={1,1,1,4,1,6,2};
184 int spin[7]={1,1,2,0,1,0,2};
185 double r0[7]={3.0,3.0,3.0,3.0,3.0,3.0,3.0};
186 double r1[7]={5.0,5.0,5.0,5.0,5.0,5.0,5.0};
187
188 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
189 setProb(value);
190 return ;
191
192}
193
194void EvtD0ToKSpi0pi0::Com_Multi(double a1[2], double a2[2], double res[2])
195{
196 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
197 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
198}
199void EvtD0ToKSpi0pi0::Com_Divide(double a1[2], double a2[2], double res[2])
200{
201 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
202 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
203 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
204}
205double EvtD0ToKSpi0pi0::CalRho4pi(double_t s)
206{
207 if(s>=1.){
208 return sqrt((s-16.*mass_Pion*mass_Pion)/s);
209 }
210 else{
211 double_t s0 = 1.2274+0.00370909/(s*s) - (0.111203)/(s) - 6.39017*s +16.8358*s*s - 21.8845*s*s*s + 11.3153*s*s*s*s;
212 double_t gam = s0*sqrt(1.0-(16.0*mass_Pion*mass_Pion));
213
214 return gam;
215 }
216}
217
218//------------base---------------------------------
219double EvtD0ToKSpi0pi0::SCADot(double a1[4], double a2[4])
220{
221 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
222 return _cal;
223}
224double EvtD0ToKSpi0pi0::barrier(int l, double sa, double sb, double sc, double r, double mass)
225{
226 double q = fabs((sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb);
227 double z;
228 z = q*r*r;
229 double sa0;
230 sa0 = mass*mass;
231 double q0 = fabs((sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb);
232 double z0 = q0*r*r;
233 double F = 0.0;
234 if(l == 0) F = 1;
235 if(l == 1) F = sqrt((1+z0)/(1+z));
236 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
237 return F;
238}
239//------------------spin-------------------------------------------
240void EvtD0ToKSpi0pi0::calt1(double daug1[4], double daug2[4], double t1[4])
241{
242 double p, pq, tmp;
243 double pa[4], qa[4];
244 for(int i=0; i<4; i++) {
245 pa[i] = daug1[i] + daug2[i];
246 qa[i] = daug1[i] - daug2[i];
247 }
248 p = SCADot(pa,pa);
249 pq = SCADot(pa,qa);
250 tmp = pq/p;
251 for(int i=0; i<4; i++) {
252 t1[i] = qa[i] - tmp*pa[i];
253 }
254}
255void EvtD0ToKSpi0pi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
256{
257 double p, r;
258 double pa[4], t1[4];
259 calt1(daug1,daug2,t1);
260 r = SCADot(t1,t1)/3.0;
261 for(int i=0; i<4; i++) {
262 pa[i] = daug1[i] + daug2[i];
263 }
264 p = SCADot(pa,pa);
265 for(int i=0; i<4; i++) {
266 for(int j=0; j<4; j++) {
267 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
268 }
269 }
270}
271//-------------------prop--------------------------------------------
272void EvtD0ToKSpi0pi0::propagator(double mass2, double mass, double width, double sx, double prop[2])
273{
274 double a[2], b[2];
275 a[0] = 1;
276 a[1] = 0;
277 b[0] = mass2-sx;
278 b[1] = -mass*width;
279 Com_Divide(a,b,prop);
280}
281double EvtD0ToKSpi0pi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
282{
283 double widm = 0.;
284 double m = sqrt(sa);
285 double tmp = sb-sc;
286 double tmp1 = sa+tmp;
287 double q = fabs(0.25*tmp1*tmp1/sa-sb);
288 double tmp2 = mass2+tmp;
289 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
290 double z = q*r2;
291 double z0 = q0*r2;
292 double t = q/q0;
293 if(l == 0) {widm = sqrt(t)*mass/m;}
294 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
295 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
296 return widm;
297}
298double EvtD0ToKSpi0pi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
299{
300 double widm = 0.;
301 double m = sqrt(sa);
302 double tmp = sb-sc;
303 double tmp1 = sa+tmp;
304 double q = fabs(0.25*tmp1*tmp1/sa-sb);
305 double tmp2 = mass2+tmp;
306 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
307 double z = q*r2;
308 double z0 = q0*r2;
309 double F = (1+z0)/(1+z);
310 double t = q/q0;
311 widm = t*sqrt(t)*mass/m*F;
312 return widm;
313}
314void EvtD0ToKSpi0pi0::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
315{
316 double a[2], b[2];
317 double mass2 = mass*mass;
318 a[0] = 1;
319 a[1] = 0;
320 b[0] = mass2-sa;
321 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
322 Com_Divide(a,b,prop);
323}
324
325void EvtD0ToKSpi0pi0::propagatorFlatte(double mass, double width, double sa, double prop[2]){
326 double q2_Pi, q2_Ka;
327 double rhoPi[2], rhoKa[2];
328
329 q2_Pi = 0.25*sa-mPi*mPi;
330 q2_Ka = 0.25*sa-mKa*mKa;
331
332 if (q2_Pi > 0) {
333 rhoPi[0] = 2.0*sqrt(q2_Pi/sa);
334 rhoPi[1] = 0.0;
335 }
336 if (q2_Pi <= 0) {
337 rhoPi[0] = 0.0;
338 rhoPi[1] = 2.0*sqrt(-q2_Pi/sa);
339 }
340
341 if (q2_Ka > 0) {
342 rhoKa[0] = 2.0*sqrt(q2_Ka/sa);
343 rhoKa[1] = 0.0;
344 }
345 if (q2_Ka <= 0) {
346 rhoKa[0] = 0.0;
347 rhoKa[1] = 2.0*sqrt(-q2_Ka/sa);
348 }
349
350 double a[2], b[2];
351 a[0] = 1;
352 a[1] = 0;
353 b[0] = mass*mass - sa + 0.165*rhoPi[1] + 0.69465*rhoKa[1];
354 b[1] = - (0.165*rhoPi[0] + 0.69465*rhoKa[0]);
355 Com_Divide(a,b,prop);
356}
357
358void EvtD0ToKSpi0pi0::rhoab(double sa, double sb, double sc, double res[2]) {
359 double tmp = sa+sb-sc;
360 double q = 0.25*tmp*tmp/sa-sb;
361 if(q>=0) {
362 res[0]=2.0*sqrt(q/sa);
363 res[1]=0.0;
364 } else {
365 res[0]=0.0;
366 res[1]=2.0*sqrt(-q/sa);
367 }
368}
369void EvtD0ToKSpi0pi0::rho4Pi(double sa, double res[2]) {
370 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
371 if(temp>=0) {
372 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
373 res[1]=0.0;
374 } else {
375 res[0]=0.0;
376 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));
377 }
378}
379
380void EvtD0ToKSpi0pi0::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
381 double f = 0.5843+1.6663*sa;
382 const double M = 0.9264;//M(f0500)
383 const double mass2 = 0.85821696; // M*M
384 const double mpi2d2 = 0.00973989245;//mass_Pion2/2 = 0.0194797849/2
385 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
386 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
387 rhoab(sa,sb,sc,rho1s);
388 rhoab(mass2,sb,sc,rho1M);
389 rho4Pi(sa,rho2s);
390 rho4Pi(mass2,rho2M);
391 Com_Divide(rho1s,rho1M,rho1);
392 Com_Divide(rho2s,rho2M,rho2);
393 double a[2], b[2];
394 a[0] = 1.0;
395 a[1] = 0.0;
396 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
397 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
398 Com_Divide(a,b,prop);
399}
400void EvtD0ToKSpi0pi0::Flatte_rhoab(double sa, double sb, double rho[2])
401{
402 double q = 1.0-(4*sb/sa);
403
404 if(q>0) {
405 rho[0]=sqrt(q);
406 rho[1]=0;
407 }
408 else if(q<0){
409 rho[0]=0;
410 rho[1]=sqrt(-q);
411 }
412}
413
414void EvtD0ToKSpi0pi0::propagator980(double mass, double sx, double *sb, double prop[2])
415{
416 double gpipi1 = 2.0/3.0*0.165;
417 double gpipi2 = 1.0/3.0*0.165;
418 double gKK1 = 0.5*0.69465;
419 double gKK2 = 0.5*0.69465;
420
421 double unit[2]={1.0};
422 double ci[2]={0,1};
423 double rho1[2];
424 Flatte_rhoab(sx,sb[0],rho1);
425 double rho2[2];
426 Flatte_rhoab(sx,sb[1],rho2);
427 double rho3[2];
428 Flatte_rhoab(sx,sb[2],rho3);
429 double rho4[2];
430 Flatte_rhoab(sx,sb[3],rho4);
431
432 double tmp1[2]={gpipi1,0};
433 double tmp11[2];
434 double tmp2[2]={gpipi2,0};
435 double tmp22[2];
436 double tmp3[2]={gKK1,0};
437 double tmp33[2];
438 double tmp4[2]={gKK2,0};
439 double tmp44[2];
440
441 Com_Multi(tmp1,rho1,tmp11);
442 Com_Multi(tmp2,rho2,tmp22);
443 Com_Multi(tmp3,rho3,tmp33);
444 Com_Multi(tmp4,rho4,tmp44);
445
446 double tmp5[2]={tmp11[0]+tmp22[0]+tmp33[0]+tmp44[0],tmp11[1]+tmp22[1]+tmp33[1]+tmp44[1]};
447 double tmp51[2];
448 Com_Multi(tmp5, ci,tmp51);
449 double tmp6[2]={mass*mass-sx-tmp51[0], -1.0*tmp51[1]};
450
451 Com_Divide( unit,tmp6, prop);
452
453}
454
455void EvtD0ToKSpi0pi0::KPiSLASS(double sa, double sb, double sc, double prop[2])
456{
457 const double m1430 = 1.441;
458 const double sa0 = 2.076481; // m1430*m1430;
459 const double w1430 = 0.193;
460 const double Lass1 = 0.25/sa0;
461 double tmp = sb-sc;
462 double tmp1 = sa0+tmp;
463 double q0 = fabs(Lass1*tmp1*tmp1-sb);
464 double tmp2 = sa+tmp;
465 double qs = 0.25*tmp2*tmp2/sa-sb;
466 double q = sqrt(qs);
467 double width = w1430*q*m1430/sqrt(sa*q0);
468 double temp_R = atan(m1430*width/(sa0-sa));
469 if(temp_R<0) temp_R += math_pi;
470 double deltaR = -109.7*math_pi/180.0 + temp_R;
471 double temp_F = atan(0.226*q/(2.0-3.8194*qs)); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
472 if(temp_F<0) temp_F += math_pi;
473 double deltaF = 0.1*math_pi/180.0 + temp_F;
474 double deltaS = deltaR + 2.0*deltaF;
475 double t1 = 0.96*sin(deltaF);
476 double t2 = sin(deltaR);
477 double CF[2], CS[2];
478 CF[0] = cos(deltaF);
479 CF[1] = sin(deltaF);
480 CS[0] = cos(deltaS);
481 CS[1] = sin(deltaS);
482 prop[0] = t1*CF[0] + t2*CS[0];
483 prop[1] = t1*CF[1] + t2*CS[1];
484}
485//------------GS---used by rho----------------------------
486void EvtD0ToKSpi0pi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
487{
488 double a[2], b[2];
489 double tmp = sb-sc;
490 double tmp1 = sa+tmp;
491 double q2 = fabs(0.25*tmp1*tmp1/sa-sb);
492
493 double tmp2 = mass2+tmp;
494 double q02 = fabs(0.25*tmp2*tmp2/mass2-sb);
495
496 double q = sqrt(q2);
497 double q0 = sqrt(q02);
498 double m = sqrt(sa);
499 double q03 = q0*q02;
500 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
501
502 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
503 double h0 = GS1*q0/mass*tmp3;
504 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
505 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
506 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
507
508 a[0] = 1.0+d*width/mass;
509 a[1] = 0.0;
510 b[0] = mass2-sa+width*f;
511 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
512 Com_Divide(a,b,prop);
513}
514
515double EvtD0ToKSpi0pi0::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
516 double pR[4], pD[4];
517 double temp_PDF, v_re;
518 temp_PDF = 0.0;
519 v_re = 0.0;
520 double B[2], s1, s2, s3, sR, sD;
521 for(int i=0; i<4; i++){
522 pR[i] = P1[i] + P2[i];
523 pD[i] = pR[i] + P3[i];
524 }
525 s1 = SCADot(P1,P1);
526 s2 = SCADot(P2,P2);
527 s3 = SCADot(P3,P3);
528 sR = SCADot(pR,pR);
529 sD = SCADot(pD,pD);
530 int G[4][4];
531 for(int i=0; i!=4; i++){
532 for(int j=0; j!=4; j++){
533 if(i==j){
534 if(i==0) G[i][j] = 1;
535 else G[i][j] = -1;
536 }
537 else G[i][j] = 0;
538 }
539 }
540 if(Ang == 0){
541 B[0] = 1;
542 B[1] = 1;
543 temp_PDF = 1;
544 }
545 if(Ang == 1){
546 B[0] = barrier(1,sR,s1,s2,3.0,mass);
547 B[1] = barrier(1,sD,sR,s3,5.0,mD0);
548 double t1[4], T1[4];
549 calt1(P1,P2,t1);
550 calt1(pR,P3,T1);
551 temp_PDF = 0;
552 for(int i=0; i<4; i++){
553 temp_PDF += t1[i]*T1[i]*G[i][i];
554 }
555 }
556 if(Ang == 2){
557 B[0] = barrier(2,sR,s1,s2,3.0,mass);
558 B[1] = barrier(2,sD,sR,s3,5.0,mD0);
559 double t2[4][4], T2[4][4];
560 calt2(P1,P2,t2);
561 calt2(pR,P3,T2);
562 temp_PDF = 0;
563 for(int i=0; i<4; i++){
564 for(int j=0; j<4; j++){
565 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
566 }
567 }
568 }
569 v_re = temp_PDF*B[0]*B[1];
570 return v_re;
571}
572
573void EvtD0ToKSpi0pi0::rhoMTX(int i, int j, double s, double Rho[2]){
574
575 double rhoijx;
576 double rhoijy;
577 double mpi = 0.13957;
578 if(i==j && i==0 ){
579 double m2 = 0.13957*0.13957;
580 if((1-(4*m2)/s)>0){
581 rhoijx = sqrt(1.0f - (4*m2)/s);
582 rhoijy = 0;
583 }else{
584 rhoijy = sqrt((4*m2)/s - 1.0f);
585 rhoijx = 0;
586 }
587 }
588 if(i==j && i==1 ){
589 double m2 = 0.49368*0.49368;
590 if((1-(4*m2)/s)>0){
591 rhoijx = sqrt(1.0f - (4*m2)/s);
592 rhoijy = 0;
593 }else{
594 rhoijy = sqrt((4*m2)/s - 1.0f);
595 rhoijx = 0;
596 }
597 }
598 if(i==j && i==2){
599 rhoijx = CalRho4pi(s);
600 rhoijy = 0;
601 }
602 if(i==j && i==3){
603 double m2 = 0.547862*0.547862;
604 if((1-(4*m2)/s)>0){
605 rhoijx = sqrt(1.0f - (4*m2)/s);
606 rhoijy = 0;
607 }else{
608 rhoijy = sqrt((4*m2)/s - 1.0f);
609 rhoijx = 0;
610 }
611 }
612 if(i==j && i==4){
613 double m_1 = 0.547862;
614 double m_2 = 0.95778;
615 double mp2 = (m_1+m_2)*(m_1+m_2);
616 double mm2 = (m_1-m_2)*(m_1-m_2);
617 if((1-mp2/s)>0){
618 rhoijx = sqrt(1.0f - mp2/s);
619 rhoijy = 0;
620 }else{
621 rhoijy = sqrt(mp2/s - 1.0f);
622 rhoijx = 0;
623 }
624 }
625
626 if(i!=j){
627 rhoijx = 0;
628 rhoijy = 0;
629 }
630 Rho[0] = rhoijx;
631 Rho[1] = rhoijy;
632
633}
634
635
636/* Kij = (sum_a gigj/(ma^2-s) + fij(1-s0)/(s-s0))((s-0.5sampi2)(1-sa0)/(s-sa0)) */
637void EvtD0ToKSpi0pi0::KMTX(int i, int j, double s,double KM[2]){
638
639 double Kijx;
640 double Kijy;
641 double mpi = 0.13957;
642 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
643
644 double g1[5] = { 0.22889,-0.55377, 0.00000,-0.39899,-0.34639};
645 double g2[5] = { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503};
646 double g3[5] = { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681};
647 double g4[5] = { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984};
648 double g5[5] = { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358};
649
650 double f1[5] = { 0.23399, 0.15044,-0.20545, 0.32825, 0.35412};
651
652 double upimag[5] = { 0,0,0,0,0};
653
654 for(int k=0; k<5; k++){
655 upimag[k] = 0;
656 }
657 double ss0 = -3.92637;
658 double sA = 1.0;//v1
659 double sA0 = -0.15;
660
661 if(i==0||j==0){
662 Kijx = (g1[i]*g1[j]/(m[0]*m[0]-s) + g2[i]*g2[j]/(m[1]*m[1]-s) + g3[i]*g3[j]/(m[2]*m[2]-s) + g4[i]*g4[j]/(m[3]*m[3]-s) + g5[i]*g5[j]/(m[4]*m[4]-s)+f1[j]*(1-ss0)/(s-ss0))*(1-sA0)/(s-sA0)*(s-sA*mpi*mpi*0.5);
663 Kijy = (g1[i]*g1[j]*upimag[0] + g2[i]*g2[j]*upimag[1] + g3[i]*g3[j]*upimag[2] + g4[i]*g4[j]*upimag[3] + g5[i]*g5[j]*upimag[4])*(1-sA0)/(s-sA0)*(s-sA*mpi*mpi*0.5);
664 }
665
666 else{
667 Kijx = (g1[i]*g1[j]/(m[0]*m[0]-s) + g2[i]*g2[j]/(m[1]*m[1]-s) + g3[i]*g3[j]/(m[2]*m[2]-s) + g4[i]*g4[j]/(m[3]*m[3]-s) + g5[i]*g5[j]/(m[4]*m[4]-s))*(1-sA0)/(s-sA0)*(s-sA*mpi*mpi*0.5);
668 Kijy = (g1[i]*g1[j]*upimag[0] + g2[i]*g2[j]*upimag[1] + g3[i]*g3[j]*upimag[2] + g4[i]*g4[j]*upimag[3] + g5[i]*g5[j]*upimag[4])*(1-sA0)/(s-sA0)*(s-sA*mpi*mpi*0.5);
669 }
670
671 KM[0] = Kijx;
672 KM[1] = Kijy;
673}
674
675
676void EvtD0ToKSpi0pi0::IMTX(int i, int j, double IMTX[2]){
677
678 double Iijx;
679 double Iijy;
680 if(i==j){
681 Iijx = 1;
682 Iijy = 0;
683 }else{
684 Iijx = 0;
685 Iijy = 0;
686 }
687 IMTX[0] = Iijx;
688 IMTX[1] = Iijy;
689
690}
691
692/* F = I - i*K*rho */
693void EvtD0ToKSpi0pi0::FMTX(double Kijx, double Kijy, double rhojjx, double rhojjy, int i, int j, double FM[2]){
694
695 double Fijx;
696 double Fijy;
697
698 double tmpx = Kijx*rhojjx - Kijy*rhojjy;
699 double tmpy = Kijy*rhojjx + Kijx*rhojjy;
700
701 double imtx[2];
702 IMTX(i,j,imtx);
703 Fijx = imtx[0]+tmpy;
704 Fijy = -tmpx;
705
706 FM[0] = Fijx;
707 FM[1] = Fijy;
708
709}
710
711void EvtD0ToKSpi0pi0::PVTR(int ID, double s, double PV[2]){
712
713 double VPix;
714 double VPiy;
715 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
716
717 double g[5][5] = {{ 0.22889,-0.55377, 0.00000,-0.39899,-0.34639},
718 { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503},
719 { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681},
720 { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984},
721 { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358}};
722
723 double betax[5], betay[5], fprodx[5], fprody[5];
724
725 betax[0] = 8.5*cos( 68.5*3.1415926/180.0); betay[0] = 8.5*sin( 68.5*3.1415926/180.0);
726 betax[1] = 12.2*cos( 24.0*3.1415926/180.0); betay[1] = 12.2*sin( 24.0*3.1415926/180.0);
727 betax[2] = 29.2*cos( -0.1*3.1415926/180.0); betay[2] = 29.2*sin( -0.1*3.1415926/180.0);
728 betax[3] = 10.8*cos(-51.9*3.1415926/180.0); betay[3] = 10.8*sin(-51.9*3.1415926/180.0);
729 betax[4] = 0.0; betay[4] = 0.0;
730
731 fprodx[0] = 8.0*cos(-126.0*3.1415926/180.0); fprody[0] = 8.0*sin(-126.0*3.1415926/180.0);
732 fprodx[1] = 26.3*cos(-152.3*3.1415926/180.0); fprody[1] = 26.3*sin(-152.3*3.1415926/180.0);
733 fprodx[2] = 33.0*cos( -93.2*3.1415926/180.0); fprody[2] = 33.0*sin( -93.2*3.1415926/180.0);
734 fprodx[3] = 26.2*cos(-121.4*3.1415926/180.0); fprody[3] = 26.2*sin(-121.4*3.1415926/180.0);
735 fprodx[4] = 0.0; fprody[4] = 0.0;
736
737 double V0x = 0.0, V0y = 0.0, V1x = 0.0, V1y = 0.0;
738 double s0_prod = -0.07;
739
740 for(int k=0;k<5;k++) {
741 V0x += betax[k]*g[k][ID]/(m[k]*m[k]-s);
742 V0y += betay[k]*g[k][ID]/(m[k]*m[k]-s);
743 }
744 V1x += (1.-s0_prod)/(s-s0_prod)*fprodx[ID];
745 V1y += (1.-s0_prod)/(s-s0_prod)*fprody[ID];
746
747 VPix = V0x+V1x;
748 VPiy = V0y+V1y;
749
750 PV[0] = VPix;
751 PV[1] = VPiy;
752
753}
754
755/* inverse for Matrix (I - i*rho*K) using LUP */
756void EvtD0ToKSpi0pi0::FINVMTX(double s, double *FINVx, double *FINVy){
757
758 int P[5] = { 0,1,2,3,4};
759
760 double Fx[5][5];
761 double Fy[5][5];
762
763 double Ux[5][5];
764 double Uy[5][5];
765 double Lx[5][5];
766 double Ly[5][5];
767
768 double UIx[5][5];
769 double UIy[5][5];
770 double LIx[5][5];
771 double LIy[5][5];
772
773 double rho[2];
774 double KM[2];
775 for(int k=0; k<5; k++){
776 rhoMTX(k,k,s,rho);
777 double rhokkx = rho[0];
778 double rhokky = rho[1];
779 Ux[k][k] = rhokkx;
780 Uy[k][k] = rhokky;
781 for(int l=k; l<5; l++){
782 KMTX(k,l,s,KM);
783 double Kklx = KM[0];
784 double Kkly = KM[1];
785 Lx[k][l] = Kklx;
786 Ly[k][l] = Kkly;
787 Lx[l][k] = Lx[k][l];
788 Ly[l][k] = Ly[k][l];
789 }
790 }
791
792 double AA[2];
793 for(int k=0; k<5; k++){
794 for(int l=0; l<5; l++){
795 FMTX(Lx[k][l],Ly[k][l],Ux[l][l],Uy[l][l],k,l,AA);
796 double Fklx = AA[0];
797 double Fkly = AA[1];
798 Fx[k][l] = Fklx;
799 Fy[k][l] = Fkly;
800 }
801 }
802
803 for(int k=0; k<5; k++){
804 double tmprM = (Fx[k][k]*Fx[k][k]+Fy[k][k]*Fy[k][k]);
805 int tmpID = 0;
806 for(int l=k; l<5; l++){
807 double tmprF = (Fx[l][k]*Fx[l][k]+Fy[l][k]*Fy[l][k]);
808 if(tmprM<=tmprF){
809 tmprM = tmprF;
810 tmpID = l;
811 }
812 }
813
814 int tmpP = P[k];
815 P[k] = P[tmpID];
816 P[tmpID] = tmpP;
817
818 for(int l=0; l<5; l++){
819
820 double tmpFx = Fx[k][l];
821 double tmpFy = Fy[k][l];
822
823 Fx[k][l] = Fx[tmpID][l];
824 Fy[k][l] = Fy[tmpID][l];
825
826 Fx[tmpID][l] = tmpFx;
827 Fy[tmpID][l] = tmpFy;
828
829 }
830 for(int l=k+1; l<5; l++){
831 double rFkk = Fx[k][k]*Fx[k][k] + Fy[k][k]*Fy[k][k];
832 double Fxlk = Fx[l][k];
833 double Fylk = Fy[l][k];
834 double Fxkk = Fx[k][k];
835 double Fykk = Fy[k][k];
836 Fx[l][k] = (Fxlk*Fxkk + Fylk*Fykk)/rFkk;
837 Fy[l][k] = (Fylk*Fxkk - Fxlk*Fykk)/rFkk;
838 for(int m=k+1; m<5; m++){
839 Fx[l][m] = Fx[l][m] - (Fx[l][k]*Fx[k][m] - Fy[l][k]*Fy[k][m]);
840 Fy[l][m] = Fy[l][m] - (Fx[l][k]*Fy[k][m] + Fy[l][k]*Fx[k][m]);
841 }
842 }
843 }
844
845 for(int k=0; k<5; k++){
846 for(int l=0; l<5 ;l++){
847 if(k==l){
848 Lx[k][k] = 1;
849 Ly[k][k] = 0;
850 Ux[k][k] = Fx[k][k];
851 Uy[k][k] = Fy[k][k];
852 }
853 if(k>l){
854 Lx[k][l] = Fx[k][l];
855 Ly[k][l] = Fy[k][l];
856 Ux[k][l] = 0;
857 Uy[k][l] = 0;
858 }
859 if(k<l){
860 Ux[k][l] = Fx[k][l];
861 Uy[k][l] = Fy[k][l];
862 Lx[k][l] = 0;
863 Ly[k][l] = 0;
864 }
865 }
866 }
867
868 for(int k=0; k<5; k++){
869
870 LIx[k][k] = 1;
871 LIy[k][k] = 0;
872
873 double rUkk = Ux[k][k]*Ux[k][k] + Uy[k][k]*Uy[k][k];
874 UIx[k][k] = Ux[k][k]/rUkk;
875 UIy[k][k] = -1.0f * Uy[k][k]/rUkk ;
876
877 for(int l=(k+1); l<5; l++){
878 LIx[k][l] = 0;
879 LIy[k][l] = 0;
880 UIx[l][k] = 0;
881 UIy[l][k] = 0;
882 }
883 for(int l=(k-1); l>=0; l--){ // U-1
884 double sx = 0;
885 double c_sx = 0;
886 double sy = 0;
887 double c_sy = 0;
888 for(int m=l+1; m<=k; m++){
889 sx = sx - c_sx;
890 double sx_tmp = sx + Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k];
891 c_sx = (sx_tmp - sx) - (Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k]);
892 sx = sx_tmp;
893
894 sy = sy - c_sy;
895 double sy_tmp = sy + Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k];
896 c_sy = (sy_tmp - sy) - (Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k]);
897 sy = sy_tmp;
898 }
899 UIx[l][k] = -1.0f * (UIx[l][l]*sx - UIy[l][l]*sy);
900 UIy[l][k] = -1.0f * (UIy[l][l]*sx + UIx[l][l]*sy);
901 }
902
903 for(int l=k+1; l<5; l++){ // L-1
904 double sx = 0;
905 double c_sx = 0;
906 double sy = 0;
907 double c_sy = 0;
908 for(int m=k; m<l; m++){
909 sx = sx - c_sx;
910 double sx_tmp = sx + Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k];
911 c_sx = (sx_tmp - sx) - (Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k]);
912 sx = sx_tmp;
913
914 sy = sy - c_sy;
915 double sy_tmp = sy + Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k];
916 c_sy = (sy_tmp - sy) - (Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k]);
917 sy = sy_tmp;
918 }
919 LIx[l][k] = -1.0f * sx;
920 LIy[l][k] = -1.0f * sy;
921 }
922 }
923 for(int m=0; m<5; m++){
924 double resX = 0;
925 double c_resX = 0;
926 double resY = 0;
927 double c_resY = 0;
928 for(int k=0; k<5; k++){
929 for(int l=0; l<5; l++){
930 double Plm = 0;
931 if(P[l] == m) Plm = 1;
932
933 resX = resX - c_resX;
934 double resX_tmp = resX + (UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm;
935 c_resX = (resX_tmp - resX) - ((UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm);
936 resX = resX_tmp;
937
938 resY = resY - c_resY;
939 double resY_tmp = resY + (UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm;
940 c_resY = (resY_tmp - resY) - ((UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm);
941 resY = resY_tmp;
942 }
943 }
944 FINVx[m] = resX;
945 FINVy[m] = resY;
946 }
947}
948
949void EvtD0ToKSpi0pi0::Fvector(double sa, double s0, double Fv[2]){
950
951 double outputx = 0;
952 double outputy = 0;
953
954 double FINVx[5] = {0,0,0,0,0};
955 double FINVy[5] = {0,0,0,0,0};
956
957 FINVMTX(sa,FINVx,FINVy);
958
959 double resx = 0;
960 double c_resx = 0;
961 double resy = 0;
962 double c_resy = 0;
963 double pv[2];
964 for(int j=0; j<5; j++){
965 PVTR(j,sa,pv);
966 double Plx = pv[0];
967 double Ply = pv[1];
968 resx = resx - c_resx;
969 double resx_tmp = resx + (FINVx[j]*Plx - FINVy[j]*Ply);
970 c_resx = (resx_tmp - resx) - (FINVx[j]*Plx - FINVy[j]*Ply);
971 resx = resx_tmp;
972
973 resy = resy - c_resy;
974 double resy_tmp = resy + (FINVx[j]*Ply + FINVy[j]*Plx);
975 c_resy = (resy_tmp - resy) - (FINVx[j]*Ply + FINVy[j]*Plx);
976 resy = resy_tmp;
977 }
978 outputx = resx;
979 outputy = resy;
980 Fv[0] = outputx;
981 Fv[1] = outputy;
982
983}
984void EvtD0ToKSpi0pi0::calEva(double* Ks0, double* Pi01, double* Pi02, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result , double*r0 ,double*r1)
985{
986
987 double P12[4], P23[4], P13[4];
988 double cof[2], amp_PDF[2], PDF[2];
989 double Amp_KPiS[2];
990 double Fv[2];
991 double s1, s2, s3;
992 double Realpipis, Imagpipis;
993
994 double s12, s13, s23;
995 for(int i=0; i<4; i++){
996 P12[i] = Ks0[i] + Pi01[i];
997 P13[i] = Ks0[i] + Pi02[i];
998 P23[i] = Pi01[i] + Pi02[i];
999 }
1000 s1 = SCADot(Ks0,Ks0);
1001 s2 = SCADot(Pi01,Pi01);
1002 s3 = SCADot(Pi02,Pi02);
1003
1004 s12 = SCADot(P12,P12);
1005 s13 = SCADot(P13,P13);
1006 s23 = SCADot(P23,P23);
1007 double pro[2], temp_PDF, pro1[2], temp_PDF1, pro2[2], temp_PDF2, amp_tmp[2];
1008 double mass1sq;
1009 amp_PDF[0] = 0;
1010 amp_PDF[1] = 0;
1011 PDF[0] = 0;
1012 PDF[1] = 0;
1013 amp_tmp[0] = 0;
1014 amp_tmp[1] = 0;
1015 for(int i=0; i<nstates; i++) {
1016 amp_tmp[0] = 0;
1017 amp_tmp[1] = 0;
1018 mass1sq = mass1[i]*mass1[i];
1019 cof[0] = amp[i]*cos(phase[i]);
1020 cof[1] = amp[i]*sin(phase[i]);
1021 temp_PDF = 0;
1022
1023 if(modetype[i] == 12){
1024 temp_PDF1 = DDalitz(Ks0, Pi01, Pi02, spin[i], mass1[i]);
1025 if(g0[i]==1) propagatorRBW(mass1[i],width1[i],s12,s1,s2,rRes2,spin[i],pro1);
1026 if(g0[i]==4) KPiSLASS(s12,s1,s2,pro1);
1027 if(g0[i]==0){
1028 pro1[0] = 1;
1029 pro1[1] = 0;
1030 }
1031
1032 temp_PDF2 = DDalitz(Ks0, Pi02, Pi01, spin[i], mass1[i]);
1033 if(g0[i]==1) propagatorRBW(mass1[i],width1[i],s13,s1,s3,rRes2,spin[i],pro2);
1034 if(g0[i]==4) KPiSLASS(s13,s1,s3,pro2);
1035 if(g0[i]==0){
1036 pro2[0] = 1;
1037 pro2[1] = 0;
1038 }
1039 amp_tmp[0] = temp_PDF1*pro1[0] + temp_PDF2*pro2[0];
1040 amp_tmp[1] = temp_PDF1*pro1[1] + temp_PDF2*pro2[1];
1041 }
1042
1043 if(modetype[i] == 23){
1044 temp_PDF = DDalitz(Pi01, Pi02, Ks0, spin[i], mass1[i]);
1045 if(g0[i]==6) {
1046 Fvector(s23,-0.07,Fv);
1047 Realpipis = Fv[0];
1048 Imagpipis = Fv[1];
1049 amp_tmp[0] = temp_PDF*Realpipis;
1050 amp_tmp[1] = temp_PDF*Imagpipis;
1051 }
1052 if(g0[i]==4) propagatorFlatte(mass1[i],width1[i],s23,pro); // Only for f0(980)
1053 if(g0[i]==2) propagatorRBW(mass1[i],width1[i],s23,s2,s3,rRes2,spin[i],pro);
1054 if(g0[i]==0){
1055 pro[0] = 1;
1056 pro[1] = 0;
1057 }
1058 if(g0[i]==500) propagatorsigma500(s23, s2, s3, pro);
1059 if(g0[i]!=6) amp_tmp[0] = temp_PDF*pro[0];
1060 if(g0[i]!=6) amp_tmp[1] = temp_PDF*pro[1];
1061 }
1062
1063 if(modetype[i] == 132){
1064 KPiSLASS(s13,s1,s3,Amp_KPiS);
1065 amp_tmp[0] = Amp_KPiS[0];
1066 amp_tmp[1] = Amp_KPiS[1];
1067 }
1068
1069 Com_Multi(amp_tmp,cof,amp_PDF);
1070 PDF[0] += amp_PDF[0];
1071 PDF[1] += amp_PDF[1];
1072
1073
1074 }
1075 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
1076 if(value <=0) value = 1e-20;
1077 Result = value;
1078 }
double sin(const BesAngle a)
Definition BesAngle.h:210
double cos(const BesAngle a)
Definition BesAngle.h:213
double P(RecMdcKalTrack *trk)
double mass
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
TFile * f1
TF1 * g1
int ID[no]
EvtComplex exp(const EvtComplex &c)
*******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
const double mpi
Definition Gam4pikp.cxx:47
XmlRpcServer s
****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
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA ** Rho(770) and Omega(782) are taken from CMD-2 F_pi fit *(hep-ex/9904027)
TTree * t
Definition binning.cxx:23
static const double radToDegrees
Definition EvtConst.hh:30
void getName(std::string &name)
void decay(EvtParticle *p)
EvtDecayBase * clone()
virtual ~EvtD0ToKSpi0pi0()
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)
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
double double * m2
Definition qcdloop1.h:75
const double b
Definition slope.cxx:9