BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDTopipi0pi0.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: EvtDTopipi0pi0.cc
11// the necessary file: EvtDTopipi0pi0.hh
12//
13// Description: D+ -> pi+ pi0 pi0 (BAM-00744)
14//
15// Modification history:
16//
17// Liaoyuan Dong Sat Dec 31 13:54:49 2022 Module created
18// Tue Jan 2 12:12:00 2024 Module updated
19//------------------------------------------------------------------------
20#include "TMatrixD.h"
21#include "TComplex.h"
23#include <fstream>
24#include <stdlib.h>
27#include "EvtGenBase/EvtPDL.hh"
32#include <string>
35#include <cmath>
36using std::endl;
37
39
40void EvtDTopipi0pi0::getName(std::string& model_name){
41 model_name="DTopipi0pi0";
42}
43
47
49 // check that there are 0 arguments
50 checkNArg(0);
51 checkNDaug(3);
56
57 phi[0] =9.32300;
58 phi[1] = 0.0;
59 phi[2] = 2.61960;
60 phi[3] = 0.0;
61
62 rho[0] = 1.26430;
63 rho[1] = 1.0;
64 rho[2] = 0.95403;
65 rho[3] = 1.0;
66
67
68 modetype[0]= 23;
69 modetype[1]= 12;
70 modetype[2]= 12;
71 modetype[3]= 23;
72
73/*
74 std::cout << "Initializing EvtDTopipi0pi0" << std::endl;
75 for (int i=0; i<4; i++) {
76 cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
77 }
78*/
79 width[0] = 0.1867;
80 width[1] = 0.1502;
81 width[2] = 0.4;
82 width[3] = 0.020;
83
84 mass[0] = 1.2755;
85 mass[1] = 0.7665;
86 mass[2] = 1.465;
87 mass[3] = 0.900;
88
89 mD = 1.86486;
90 mDs = 1.9683;
91 rRes = 9.0;
92 rD = 5.0;
93 metap = 0.95778;
94 mkstr = 0.89594;
95 mk0 = 0.497614;
96 mass_Kaon = 0.49368;
97 mass_Pion = 0.13957;
98 mass_Pi0 = 0.1349766;
99 math_pi = 3.1415926;
100 ma0 = 0.99;
101 Ga0 = 0.0756;
102 meta = 0.547862;
103 mass_Kaon = 0.49368;
104 mass_Eta = 0.547862;
105 mass_Etap = 0.95778;
106
107
108 GS1 = 0.636619783;
109 GS2 = 0.01860182466;
110 GS3 = 0.1591549458; // 1/(2*math_2pi)
111 GS4 = 0.00620060822; // mass_Pion2/math_pi
112
113 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
114 for (int i=0; i<4; i++) {
115 for (int j=0; j<4; j++) {
116 G[i][j] = GG[i][j];
117 }
118 }
119}
120
122 setProbMax(528.0);
123}
124
126
128 EvtVector4R D1 = p->getDaug(0)->getP4();
129 EvtVector4R D2 = p->getDaug(1)->getP4();
130 EvtVector4R D3 = p->getDaug(2)->getP4();
131
132 double P1[4], P2[4], P3[4];
133 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
134 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
135 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
136
137 double value;
138 int g0[4]={1,2,2,6};
139 int spin[4]={2,1,1,0};
140 int nstates=4;
141 double r0[4]={3,3,3,3};
142 double r1[4]={5,5,5,5};
143
144 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
145
146 setProb(value);
147 return ;
148
149}
150
151void EvtDTopipi0pi0::Com_Multi(double a1[2], double a2[2], double res[2])
152{
153 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
154 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
155}
156void EvtDTopipi0pi0::Com_Divide(double a1[2], double a2[2], double res[2])
157{
158 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
159 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
160 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
161}
162
163double EvtDTopipi0pi0::CalRho4pi(double_t s)
164{
165 if(s>=1.){
166 return sqrt((s-16.*mass_Pion*mass_Pion)/s);
167 }
168 else{
169 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;
170 double_t gam = s0*sqrt(1.0-(16.0*mass_Pion*mass_Pion));
171
172 return gam;
173 }
174}
175
176//------------base---------------------------------
177double EvtDTopipi0pi0::SCADot(double a1[4], double a2[4])
178{
179 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
180 return _cal;
181}
182double EvtDTopipi0pi0::barrier(int l, double sa, double sb, double sc, double r, double mass)
183{
184 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
185 if(q < 0) q = 1e-16;
186 double z;
187 z = q*r*r;
188 double sa0;
189 sa0 = mass*mass;
190 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
191 if(q0 < 0) q0 = 1e-16;
192 double z0 = q0*r*r;
193 double F = 0.0;
194 if(l == 0) F = 1;
195 if(l == 1) F = sqrt((1+z0)/(1+z));
196 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
197 return F;
198}
199//------------------spin-------------------------------------------
200void EvtDTopipi0pi0::calt1(double daug1[4], double daug2[4], double t1[4])
201{
202 double p, pq, tmp;
203 double pa[4], qa[4];
204 for(int i=0; i<4; i++) {
205 pa[i] = daug1[i] + daug2[i];
206 qa[i] = daug1[i] - daug2[i];
207 }
208 p = SCADot(pa,pa);
209 pq = SCADot(pa,qa);
210 tmp = pq/p;
211 for(int i=0; i<4; i++) {
212 t1[i] = qa[i] - tmp*pa[i];
213 }
214}
215void EvtDTopipi0pi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
216{
217 double p, r;
218 double pa[4], t1[4];
219 calt1(daug1,daug2,t1);
220 r = SCADot(t1,t1)/3.0;
221 for(int i=0; i<4; i++) {
222 pa[i] = daug1[i] + daug2[i];
223 }
224 p = SCADot(pa,pa);
225 for(int i=0; i<4; i++) {
226 for(int j=0; j<4; j++) {
227 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
228 }
229 }
230}
231//-------------------prop--------------------------------------------
232void EvtDTopipi0pi0::propagator(double mass2, double mass, double width, double sx, double prop[2])
233{
234 double a[2], b[2];
235 a[0] = 1;
236 a[1] = 0;
237 b[0] = mass2-sx;
238 b[1] = -mass*width;
239 Com_Divide(a,b,prop);
240}
241double EvtDTopipi0pi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
242{
243 double widm = 0.;
244 double m = sqrt(sa);
245 double tmp = sb-sc;
246 double tmp1 = sa+tmp;
247 double q = 0.25*tmp1*tmp1/sa-sb;
248 if(q<0) q = 1e-16;
249 double tmp2 = mass2+tmp;
250 double q0 = 0.25*tmp2*tmp2/mass2-sb;
251 if(q0<0) q0 = 1e-16;
252 double z = q*r2;
253 double z0 = q0*r2;
254 double t = q/q0;
255 if(l == 0) {widm = sqrt(t)*mass/m;}
256 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
257 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
258 return widm;
259}
260double EvtDTopipi0pi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
261{
262 double widm = 0.;
263 double m = sqrt(sa);
264 double tmp = sb-sc;
265 double tmp1 = sa+tmp;
266 double q = 0.25*tmp1*tmp1/sa-sb;
267 if(q<0) q = 1e-16;
268 double tmp2 = mass2+tmp;
269 double q0 = 0.25*tmp2*tmp2/mass2-sb;
270 if(q0<0) q0 = 1e-16;
271 double z = q*r2;
272 double z0 = q0*r2;
273 double F = (1+z0)/(1+z);
274 double t = q/q0;
275 widm = t*sqrt(t)*mass/m*F;
276 return widm;
277}
278void EvtDTopipi0pi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
279{
280
281 double a[2], b[2];
282 a[0] = 1;
283 a[1] = 0;
284 b[0] = mass2-sa;
285 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
286 Com_Divide(a,b,prop);
287}
288
289void EvtDTopipi0pi0::propagatorFlatte(double mass, double width, double sa, double sb, double sc, int r, double prop[2]){
290 double q, qKsK;
291 double rhoab[2], rhoKsK[2];
292 q = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
293 if(r == 0) qKsK = 0.25*sa - 0.49368*0.49368;
294 if(r == 1) qKsK = 0.25*(sa + 3.899750596e-03)*(sa + 3.899750596e-03)/sa - 0.497614*0.497614;
295 if(q>0){
296 rhoab[0] = 2*sqrt(q/sa);
297 rhoab[1] = 0;
298 }
299 if(q<0){
300 rhoab[0] = 0;
301 rhoab[1] = 2*sqrt(-q/sa);
302 }
303 if(qKsK>0){
304 rhoKsK[0] = 2*sqrt(qKsK/sa);
305 rhoKsK[1] = 0;
306 }
307 if(qKsK<0){
308 rhoKsK[0] = 0;
309 rhoKsK[1] = 2*sqrt(-qKsK/sa);
310 }
311 double a[2], b[2];
312 a[0] = 1;
313 a[1] = 0;
314 b[0] = mass*mass - sa + 0.341*rhoab[1] + 0.892*0.341*rhoKsK[1];
315 b[1] = - (0.341*rhoab[0] + 0.892*0.341*rhoKsK[0]);
316 Com_Divide(a,b,prop);
317}
318
319void EvtDTopipi0pi0::rhoab(double sa, double sb, double sc, double res[2]) {
320 double tmp = sa+sb-sc;
321 double q = 0.25*tmp*tmp/sa-sb;
322 if(q>=0) {
323 res[0]=2.0*sqrt(q/sa);
324 res[1]=0.0;
325 } else {
326 res[0]=0.0;
327 res[1]=2.0*sqrt(-q/sa);
328 }
329}
330void EvtDTopipi0pi0::rho4Pi(double sa, double res[2]) {
331 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
332 if(temp>=0) {
333 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
334 res[1]=0.0;
335 } else {
336 res[0]=0.0;
337 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));//?????????????????????
338 }
339}
340
341
342void EvtDTopipi0pi0::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
343 double f = 0.5843+1.6663*sa;
344 const double M = 0.9264;//M(f0500)
345 const double mass2 = 0.85821696; // M*M
346 const double mpi2d2 = 0.00973989245;//mass_Pion2/2 = 0.0194797849/2
347 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
348 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
349 rhoab(sa,sb,sc,rho1s);
350 rhoab(mass2,sb,sc,rho1M);
351 rho4Pi(sa,rho2s);
352 rho4Pi(mass2,rho2M);
353 Com_Divide(rho1s,rho1M,rho1);
354 Com_Divide(rho2s,rho2M,rho2);
355 double a[2], b[2];
356 a[0] = 1.0;
357 a[1] = 0.0;
358 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
359 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
360 Com_Divide(a,b,prop);
361}
362void EvtDTopipi0pi0::Flatte_rhoab(double sa, double sb, double rho[2])
363{
364 double q = 1.0-(4*sb/sa);
365
366 if(q>0) {
367 rho[0]=sqrt(q);
368 rho[1]=0;
369 }
370 else if(q<0){
371 rho[0]=0;
372 rho[1]=sqrt(-q);
373 }
374}
375
376void EvtDTopipi0pi0::propagator980(double mass, double sx, double *sb, double prop[2])
377{
378
379 double gpipi1 = 2.0/3.0*0.165;
380 double gpipi2 = 1.0/3.0*0.165;
381 double gKK1 = 0.5*0.69465;
382 double gKK2 = 0.5*0.69465;
383
384 double unit[2]={1.0};
385 double ci[2]={0,1};
386 double rho1[2];
387 Flatte_rhoab(sx,sb[0],rho1);
388 double rho2[2];
389 Flatte_rhoab(sx,sb[1],rho2);
390 double rho3[2];
391 Flatte_rhoab(sx,sb[2],rho3);
392 double rho4[2];
393 Flatte_rhoab(sx,sb[3],rho4);
394
395 double tmp1[2]={gpipi1,0};
396 double tmp11[2];
397 double tmp2[2]={gpipi2,0};
398 double tmp22[2];
399 double tmp3[2]={gKK1,0};
400 double tmp33[2];
401 double tmp4[2]={gKK2,0};
402 double tmp44[2];
403
404 Com_Multi(tmp1,rho1,tmp11);
405 Com_Multi(tmp2,rho2,tmp22);
406 Com_Multi(tmp3,rho3,tmp33);
407 Com_Multi(tmp4,rho4,tmp44);
408
409 double tmp5[2]={tmp11[0]+tmp22[0]+tmp33[0]+tmp44[0],tmp11[1]+tmp22[1]+tmp33[1]+tmp44[1]};
410 double tmp51[2];
411 Com_Multi(tmp5, ci,tmp51);
412 double tmp6[2]={mass*mass-sx-tmp51[0], -1.0*tmp51[1]};
413
414 Com_Divide( unit,tmp6, prop);
415
416}
417
418//------------GS---used by rho----------------------------
419void EvtDTopipi0pi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
420{
421 double a[2], b[2];
422 double tmp = sb-sc;
423 double tmp1 = sa+tmp;
424 double q2 = 0.25*tmp1*tmp1/sa-sb;
425 if(q2<0) q2 = 1e-16;
426
427 double tmp2 = mass2+tmp;
428 double q02 = 0.25*tmp2*tmp2/mass2-sb;
429 if(q02<0) q02 = 1e-16;
430
431 double q = sqrt(q2);
432 double q0 = sqrt(q02);
433 double m = sqrt(sa);
434 double q03 = q0*q02;
435 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
436
437 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
438 double h0 = GS1*q0/mass*tmp3;
439 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
440 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
441 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
442
443 a[0] = 1.0+d*width/mass;
444 a[1] = 0.0;
445 b[0] = mass2-sa+width*f;
446 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
447 Com_Divide(a,b,prop);
448}
449
450double EvtDTopipi0pi0::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass, double rRES0, double rRES1){
451 double pR[4], pD[4];
452 double temp_PDF, v_re;
453 temp_PDF = 0.0;
454 v_re = 0.0;
455 double B[2], s1, s2, s3, sR, sD;
456 for(int i=0; i<4; i++){
457 pR[i] = P1[i] + P2[i];
458 pD[i] = pR[i] + P3[i];
459 }
460 s1 = SCADot(P1,P1);
461 s2 = SCADot(P2,P2);
462 s3 = SCADot(P3,P3);
463 sR = SCADot(pR,pR);
464 sD = SCADot(pD,pD);
465 int G[4][4];
466 for(int i=0; i!=4; i++){
467 for(int j=0; j!=4; j++){
468 if(i==j){
469 if(i==0) G[i][j] = 1;
470 else G[i][j] = -1;
471 }
472 else G[i][j] = 0;
473 }
474 }
475 if(Ang == 0){
476 B[0] = 1;
477 B[1] = 1;
478 temp_PDF = 1;
479 }
480 if(Ang == 1){
481 B[0] = barrier(1,sR,s1,s2,rRES0,mass);
482 B[1] = barrier(1,sD,sR,s3,rRES1,1.86966);
483 double t1[4], T1[4];
484 calt1(P1,P2,t1);
485 calt1(pR,P3,T1);
486 temp_PDF = 0;
487 for(int i=0; i<4; i++){
488 temp_PDF += t1[i]*T1[i]*G[i][i];
489 }
490 }
491 if(Ang == 2){
492 B[0] = barrier(2,sR,s1,s2,rRES0,mass);
493 B[1] = barrier(2,sD,sR,s3,rRES1,1.86966);
494 double t2[4][4], T2[4][4];
495 calt2(P1,P2,t2);
496 calt2(pR,P3,T2);
497 temp_PDF = 0;
498 for(int i=0; i<4; i++){
499 for(int j=0; j<4; j++){
500 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
501 }
502 }
503 }
504 v_re = temp_PDF*B[0]*B[1];
505 return v_re;
506}
507
508TComplex EvtDTopipi0pi0::ResonanceSkm(Double_t& m2)
509{
510 Double_t g11 = 0.22889, g12 = -0.55377, g13 = 0.00, g14 = -0.39899, g15 = -0.34639;
511 Double_t g21 = 0.94128, g22 = 0.55095, g23 = 0.00, g24 = 0.39065, g25 = 0.31503;
512 Double_t g31 = 0.36856, g32 = 0.23888, g33 = 0.55639, g34 = 0.18340, g35 = 0.18681;
513 Double_t g41 = 0.33650, g42 = 0.40907, g43 = 0.85679, g44 = 0.19906, g45 = -0.00984;
514 Double_t g51 = 0.18171, g52 = -0.17558, g53 = -0.79658, g54 = -0.00355, g55 = 0.22358;
515
516 Double_t fs11 = 0.23399, fs12 = 0.15044, fs13 =-0.20545, fs14 = 0.32825, fs15 = 0.35412;
517 Double_t m_1 = 0.65100, m_2 = 1.20360, m_3 = 1.55817, m_4 = 1.21000, m_5 = 1.82206;
518 Double_t ss0 = -3.92637;
519 Double_t sp0 = -0.07;//v1
520 double_t sA0 = -0.15;
521 double_t sA = 1.0;
522
523 double_t km11 = (g11*g11/(m_1*m_1-m2)+g21*g21/(m_2*m_2-m2)+g31*g31/(m_3*m_3-m2)+g41*g41/(m_4*m_4-m2)+g51*g51/(m_5*m_5-m2)+fs11*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
524 double_t km12 = (g11*g12/(m_1*m_1-m2)+g21*g22/(m_2*m_2-m2)+g31*g32/(m_3*m_3-m2)+g41*g42/(m_4*m_4-m2)+g51*g52/(m_5*m_5-m2)+fs12*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
525 double_t km13 = (g11*g13/(m_1*m_1-m2)+g21*g23/(m_2*m_2-m2)+g31*g33/(m_3*m_3-m2)+g41*g43/(m_4*m_4-m2)+g51*g53/(m_5*m_5-m2)+fs13*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
526 double_t km14 = (g11*g14/(m_1*m_1-m2)+g21*g24/(m_2*m_2-m2)+g31*g34/(m_3*m_3-m2)+g41*g44/(m_4*m_4-m2)+g51*g54/(m_5*m_5-m2)+fs14*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
527 double_t km15 = (g11*g15/(m_1*m_1-m2)+g21*g25/(m_2*m_2-m2)+g31*g35/(m_3*m_3-m2)+g41*g45/(m_4*m_4-m2)+g51*g55/(m_5*m_5-m2)+fs15*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
528
529 double_t km21 = km12, km31 = km13, km41 = km14, km51 = km15;
530 double_t km23 = (g12*g13/(m_1*m_1-m2)+g22*g23/(m_2*m_2-m2)+g32*g33/(m_3*m_3-m2)+g42*g43/(m_4*m_4-m2)+g52*g53/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
531 double_t km24 = (g12*g14/(m_1*m_1-m2)+g22*g24/(m_2*m_2-m2)+g32*g34/(m_3*m_3-m2)+g42*g44/(m_4*m_4-m2)+g52*g54/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
532 double_t km25 = (g12*g15/(m_1*m_1-m2)+g22*g25/(m_2*m_2-m2)+g32*g35/(m_3*m_3-m2)+g42*g45/(m_4*m_4-m2)+g52*g55/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
533 double_t km32 = km23, km42 = km24, km52 = km25;
534 double_t km34 = (g13*g14/(m_1*m_1-m2)+g23*g24/(m_2*m_2-m2)+g33*g34/(m_3*m_3-m2)+g43*g44/(m_4*m_4-m2)+g53*g54/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
535 double_t km35 = (g13*g15/(m_1*m_1-m2)+g23*g25/(m_2*m_2-m2)+g33*g35/(m_3*m_3-m2)+g43*g45/(m_4*m_4-m2)+g53*g55/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
536 double_t km43 = km34, km53 = km35;
537 double_t km45 = (g14*g15/(m_1*m_1-m2)+g24*g25/(m_2*m_2-m2)+g34*g35/(m_3*m_3-m2)+g44*g45/(m_4*m_4-m2)+g54*g55/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
538 double_t km54 = km45;
539 double_t km22 = (g12*g12/(m_1*m_1-m2)+g22*g22/(m_2*m_2-m2)+g32*g32/(m_3*m_3-m2)+g42*g42/(m_4*m_4-m2)+g52*g52/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
540 double_t km33 = (g13*g13/(m_1*m_1-m2)+g23*g23/(m_2*m_2-m2)+g33*g33/(m_3*m_3-m2)+g43*g43/(m_4*m_4-m2)+g53*g53/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
541 double_t km44 = (g14*g14/(m_1*m_1-m2)+g24*g24/(m_2*m_2-m2)+g34*g34/(m_3*m_3-m2)+g44*g44/(m_4*m_4-m2)+g54*g54/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
542 double_t km55 = (g15*g15/(m_1*m_1-m2)+g25*g25/(m_2*m_2-m2)+g35*g35/(m_3*m_3-m2)+g45*g45/(m_4*m_4-m2)+g55*g55/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
543
544 TComplex fp11(1.41620*cos(5.90540), 1.41620*sin(5.90540));
545 TComplex fp12(15.3640*cos(2.22510), 15.3640*sin(2.22510));
546 TComplex fp13(8.2297*cos(2.30550), 8.2297*sin(2.30550));
547 TComplex fp14(0.0, 0.0);
548 TComplex fp15(0.0, 0.0);
549
550 TComplex beta1(3.70140*cos(5.0818), 3.70140*sin(5.0818));
551 TComplex beta2(4.04470*cos(0.13605), 4.04470*sin(0.13605));
552 TComplex beta3(3.4591*cos(-5.6845), 3.4591*sin(-5.6845));
553 TComplex beta4(0.0, 0.0);
554 TComplex beta5(0.0, 0.0);
555
556
557
558
559 TComplex P1 = fp11*(1-sp0)/(m2-sp0)+beta1*g11/(m_1*m_1-m2)+beta2*g21/(m_2*m_2-m2)+beta3*g31/(m_3*m_3-m2)+beta4*g41/(m_4*m_4-m2)+beta5*g51/(m_5*m_5-m2);
560 TComplex P2 = fp12*(1-sp0)/(m2-sp0)+beta1*g12/(m_1*m_1-m2)+beta2*g22/(m_2*m_2-m2)+beta3*g32/(m_3*m_3-m2)+beta4*g42/(m_4*m_4-m2)+beta5*g52/(m_5*m_5-m2);
561 TComplex P3 = fp13*(1-sp0)/(m2-sp0)+beta1*g13/(m_1*m_1-m2)+beta2*g23/(m_2*m_2-m2)+beta3*g33/(m_3*m_3-m2)+beta4*g43/(m_4*m_4-m2)+beta5*g53/(m_5*m_5-m2);
562 TComplex P4 = fp14*(1-sp0)/(m2-sp0)+beta1*g14/(m_1*m_1-m2)+beta2*g24/(m_2*m_2-m2)+beta3*g34/(m_3*m_3-m2)+beta4*g44/(m_4*m_4-m2)+beta5*g54/(m_5*m_5-m2);
563 TComplex P5 = fp15*(1-sp0)/(m2-sp0)+beta1*g15/(m_1*m_1-m2)+beta2*g25/(m_2*m_2-m2)+beta3*g35/(m_3*m_3-m2)+beta4*g45/(m_4*m_4-m2)+beta5*g55/(m_5*m_5-m2);
564
565 TMatrixD MI(5,5), MA(5,5), MA_invt(5,5), MB(5,5), KM(5,5), RhoA(5,5), RhoB(5,5), MRe(5,5), MIm(5,5);
566 MI.UnitMatrix();
567 RhoA.UnitMatrix();
568 RhoB.UnitMatrix();
569 TMatrixDRow(KM,0)(0) = km11;TMatrixDRow(KM,1)(0) = km21;TMatrixDRow(KM,2)(0) = km31;TMatrixDRow(KM,3)(0) = km41;TMatrixDRow(KM,4)(0)= km51;
570 TMatrixDRow(KM,0)(1) = km12;TMatrixDRow(KM,1)(1) = km22;TMatrixDRow(KM,2)(1) = km32;TMatrixDRow(KM,3)(1) = km42;TMatrixDRow(KM,4)(1)= km52;
571 TMatrixDRow(KM,0)(2) = km13;TMatrixDRow(KM,1)(2) = km23;TMatrixDRow(KM,2)(2) = km33;TMatrixDRow(KM,3)(2) = km43;TMatrixDRow(KM,4)(2)= km53;
572 TMatrixDRow(KM,0)(3) = km14;TMatrixDRow(KM,1)(3) = km24;TMatrixDRow(KM,2)(3) = km34;TMatrixDRow(KM,3)(3) = km44;TMatrixDRow(KM,4)(3)= km54;
573 TMatrixDRow(KM,0)(4) = km15;TMatrixDRow(KM,1)(4) = km25;TMatrixDRow(KM,2)(4) = km35;TMatrixDRow(KM,3)(4) = km45;TMatrixDRow(KM,4)(4)= km55;
574
575 if ( (1.0-4.0*mass_Pion*mass_Pion/m2) > 0) {
576 TMatrixDRow(RhoA,0)(0) = sqrt(1.0-4.0*mass_Pion*mass_Pion/m2);
577 TMatrixDRow(RhoB,0)(0) = 0.0;
578 }
579 else {
580 TMatrixDRow(RhoA,0)(0) = 0.0;
581 TMatrixDRow(RhoB,0)(0) = sqrt(4.0*mass_Pion*mass_Pion/m2-1.0);
582 }
583
584 if ( (1.0-4.0*mass_Kaon*mass_Kaon/m2) > 0) {
585 TMatrixDRow(RhoA,1)(1) = sqrt(1.0-4.0*mass_Kaon*mass_Kaon/m2);
586 TMatrixDRow(RhoB,1)(1) = 0.0;
587 }
588 else {
589 TMatrixDRow(RhoA,1)(1) = 0.0;
590 TMatrixDRow(RhoB,1)(1) = sqrt(4.0*mass_Kaon*mass_Kaon/m2-1.0);
591 }
592
593 TMatrixDRow(RhoA,2)(2) = CalRho4pi(m2);
594 TMatrixDRow(RhoB,2)(2) = 0.0;
595 if ( (1.0-4.0*mass_Eta*mass_Eta/m2) > 0) {
596 TMatrixDRow(RhoA,3)(3) = sqrt(1.0-4.0*mass_Eta*mass_Eta/m2);
597 TMatrixDRow(RhoB,3)(3) = 0.0;
598 }else {
599 TMatrixDRow(RhoA,3)(3) = 0.0;
600 TMatrixDRow(RhoB,3)(3) = sqrt(4.0*mass_Eta*mass_Eta/m2-1.0);
601 }
602
603
604 if ( (1.0-(mass_Eta + mass_Etap)*(mass_Eta + mass_Etap)/m2) > 0) {
605 TMatrixDRow(RhoA,4)(4) = sqrt(1.0-(mass_Eta + mass_Etap)*(mass_Eta + mass_Etap)/m2);
606 TMatrixDRow(RhoB,4)(4) = 0.0;
607 }else{
608 TMatrixDRow(RhoA,4)(4) = 0.0;
609 TMatrixDRow(RhoB,4)(4) = sqrt((mass_Eta + mass_Etap)*(mass_Eta + mass_Etap)/m2-1.0);
610 }
611
612 MA = MI + KM*RhoB;
613 MB = -1.0*KM*RhoA;
614
615 MA_invt = MA;
616 MA_invt.Invert();
617
618 MRe = MA+MB*MA_invt*MB;
619 MRe.Invert();
620 MIm = MA_invt*MB*MRe;
621
622 TComplex ciR(1.0, 0.0);
623 TComplex ciM(0.0, 1.0);
624 TComplex amp;
625
626 amp = (ciR*TMatrixDRow(MRe,0)(0)-ciM*TMatrixDRow(MIm,0)(0))*P1+
627 (ciR*TMatrixDRow(MRe,0)(1)-ciM*TMatrixDRow(MIm,0)(1))*P2+
628 (ciR*TMatrixDRow(MRe,0)(2)-ciM*TMatrixDRow(MIm,0)(2))*P3+
629 (ciR*TMatrixDRow(MRe,0)(3)-ciM*TMatrixDRow(MIm,0)(3))*P4+
630 (ciR*TMatrixDRow(MRe,0)(4)-ciM*TMatrixDRow(MIm,0)(4))*P5;
631
632 return amp;
633}
634
635void EvtDTopipi0pi0::calEva(double* Pic, 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)
636{
637 double P12[4], P23[4], P13[4];
638 //double Pic[4], Pi01[4], Pi02[4], P12[4], P23[4], P13[4];
639 double cof[2], amp_PDF[2], PDF[2];
640 double spi01, spi02, scpi;
641 double Fv[2];
642 double s12, s13, s23;
643 TComplex tmpswave;
644 double Realpipis, Imagpipis;
645 for(int i=0; i<4; i++){
646 P23[i] = Pi01[i] + Pi02[i];
647 P12[i] = Pic[i] + Pi01[i];
648 P13[i] = Pic[i] + Pi02[i];
649 }
650 scpi = SCADot(Pic,Pic);
651 spi01 = SCADot(Pi01,Pi01);
652 spi02 = SCADot(Pi02,Pi02);
653 s23 = SCADot(P23,P23);
654 s12 = SCADot(P12,P12);
655 s13 = SCADot(P13,P13); //pipi0
656
657 double spi012[4]={mass_Pion*mass_Pion,spi02,mass_Kaon*mass_Kaon,mk0*mk0};
658
659
660
661 double pro[2], temp_PDF, pro1[2], temp_PDF1, pro2[2], temp_PDF2, amp_tmp[2];
662 double mass1sq;
663 amp_PDF[0] = 0;
664 amp_PDF[1] = 0;
665 PDF[0] = 0;
666 PDF[1] = 0;
667 amp_tmp[0] = 0;
668 amp_tmp[1] = 0;
669 for(int i=0; i<nstates; i++) {
670 amp_tmp[0] = 0;
671 amp_tmp[1] = 0;
672 mass1sq = mass1[i]*mass1[i];
673 cof[0] = amp[i]*cos(phase[i]);
674 cof[1] = amp[i]*sin(phase[i]);
675 temp_PDF = 0;
676 if(modetype[i] == 23){//pi0pi0
677 temp_PDF = DDalitz(Pi01, Pi02, Pic, spin[i], mass1[i], r0[i],r1[i]);
678 if(g0[i]==2) propagatorsigma500(s23,spi01,spi02,pro);
679 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s23,spi01,spi02,r0[i]*r0[i],spin[i],pro);
680 if(g0[i]==3) propagator980(mass1[i],s23,spi012,pro);
681 if(g0[i]==6) {
682 tmpswave = ResonanceSkm(s23);
683 Realpipis = tmpswave.Re();
684 Imagpipis = tmpswave.Im();
685 amp_tmp[0] = temp_PDF*Realpipis;
686 amp_tmp[1] = temp_PDF*Imagpipis;
687 }
688 if(g0[i]==0){
689 pro[0] = 1;
690 pro[1] = 0;
691 }
692 if(g0[i]!=6) amp_tmp[0] = temp_PDF*pro[0];
693 if(g0[i]!=6) amp_tmp[1] = temp_PDF*pro[1];
694 }
695 if(modetype[i] == 12){//pipi0
696 temp_PDF1 = DDalitz(Pic, Pi01, Pi02, spin[i], mass1[i], r0[i], r1[i]);
697 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,scpi,spi01,r0[i]*r0[i],spin[i],pro1);
698 if(g0[i]==2) propagatorGS(mass1sq,mass1[i],width1[i],s12,scpi,spi01,r0[i]*r0[i],pro1);
699 if(g0[i]==0){
700 pro1[0] = 1;
701 pro1[1] = 0;
702 }
703
704 temp_PDF2 = DDalitz(Pic, Pi02, Pi01, spin[i], mass1[i], r0[i], r1[i]);
705 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,scpi,spi02,r0[i]*r0[i],spin[i],pro2);
706 if(g0[i]==2) propagatorGS(mass1sq,mass1[i],width1[i],s13,scpi,spi02,r0[i]*r0[i],pro2);
707 if(g0[i]==0){
708 pro2[0] = 1;
709 pro2[1] = 0;
710 }
711 amp_tmp[0] = temp_PDF1*pro1[0] + temp_PDF2*pro2[0];
712
713 amp_tmp[1] = temp_PDF1*pro1[1] + temp_PDF2*pro2[1];
714 }
715 Com_Multi(amp_tmp,cof,amp_PDF);
716 PDF[0] += amp_PDF[0];
717 PDF[1] += amp_PDF[1];
718 }
719 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
720 if(value <=0) value = 1e-20;
721 Result = value;
722}
723
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")
TF1 * g1
EvtComplex exp(const EvtComplex &c)
const double mass_Pion
*******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
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
TTree * t
Definition binning.cxx:23
virtual ~EvtDTopipi0pi0()
void decay(EvtParticle *p)
EvtDecayBase * clone()
void getName(std::string &name)
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