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