BOSS 7.1.0
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//
12// Description: D+ -> pi+ pi0 pi0
13//
14// Modification history:
15//
16// Liaoyuan Dong Feb 17 02:31:49 2023 Module created
17//
18//------------------------------------------------------------------------
19#include "EvtGenBase/EvtPatches.hh"
20#include <fstream>
21#include <stdlib.h>
22#include "EvtGenBase/EvtParticle.hh"
23#include "EvtGenBase/EvtGenKine.hh"
24#include "EvtGenBase/EvtPDL.hh"
25#include "EvtGenBase/EvtReport.hh"
26#include "EvtGenBase/EvtResonance.hh"
27#include "EvtGenBase/EvtResonance2.hh"
28#include "EvtGenModels/EvtDTopipi0pi0.hh"
29#include <string>
30#include "EvtGenBase/EvtConst.hh"
31#include "EvtGenBase/EvtDecayTable.hh"
32using std::endl;
33
34EvtDTopipi0pi0::~EvtDTopipi0pi0() {}
35
36void EvtDTopipi0pi0::getName(std::string& model_name){
37 model_name="DTopipi0pi0";
38}
39
40EvtDecayBase* EvtDTopipi0pi0::clone(){
41 return new EvtDTopipi0pi0;
42}
43
44void EvtDTopipi0pi0::init(){
45 // check that there are 0 arguments
46 checkNArg(0);
47 checkNDaug(3);
48 checkSpinParent(EvtSpinType::SCALAR);
49 checkSpinDaughter(0,EvtSpinType::SCALAR);
50 checkSpinDaughter(1,EvtSpinType::SCALAR);
51 checkSpinDaughter(2,EvtSpinType::SCALAR);
52
53 phi[0] = 0.286656178602813;
54 phi[1] = -0.888625014448000;
55 phi[2] = -2.453879234281807;
56 phi[3] = -2.975162215957068;
57 phi[4] = 0.0;
58 phi[5] = 2.704379334027896;
59
60 rho[0] = 3.492106794347418;
61 rho[1] = 0.727208768061443;
62 rho[2] = 0.419996095917039;
63 rho[3] = 1.057515286243158;
64 rho[4] = 1.0;
65 rho[5] = 0.912113205640503;
66
67 modetype[0]= 23;
68 modetype[1]= 23;
69 modetype[2]= 23;
70 modetype[3]= 23;
71 modetype[4]= 12;
72 modetype[5]= 12;
73
74 //std::cout << "Initializing EvtDTopipi0pi0" << std::endl;
75 //for (int i=0; i<6; i++) {
76 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
77 //}
78
79 width[0] = 0.534;
80 width[1] = 0.0674;
81 width[2] = 0.265;
82 width[3] = 0.1867;
83 width[4] = 0.1502;
84 width[5] = 0.4;
85
86 mass[0] = 0.526;
87 mass[1] = 0.965;
88 mass[2] = 1.35;
89 mass[3] = 1.2755;
90 mass[4] = 0.7665;
91 mass[5] = 1.465;
92
93 mD = 1.86486;
94 mDs = 1.9683;
95 rRes = 9.0;
96 rD = 5.0;
97 metap = 0.95778;
98 mkstr = 0.89594;
99 mk0 = 0.497614;
100 mass_Kaon = 0.49368;
101 mass_Pion = 0.13957;
102 mass_Pi0 = 0.1349766;
103 math_pi = 3.1415926;
104 ma0 = 0.99;
105 Ga0 = 0.0756;
106 meta = 0.547862;
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
121void EvtDTopipi0pi0::initProbMax() {
122 setProbMax(495);//MAXprob = 494.1923984676
123}
124
125void EvtDTopipi0pi0::decay( EvtParticle *p){
126/*
127 double maxprob = 0.0;
128 for(int ir=0;ir<=60000000;ir++){
129 p->initializePhaseSpace(getNDaug(),getDaugs());
130 EvtVector4R D1 = p->getDaug(0)->getP4();
131 EvtVector4R D2 = p->getDaug(1)->getP4();
132 EvtVector4R D3 = p->getDaug(2)->getP4();
133
134 double P1[4], P2[4], P3[4];
135 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
136 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
137 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
138
139 double value;
140 int g0[6]={2,3,1,1,2,2};
141 int spin[6]={0,0,0,2,1,1};
142 int nstates=6;
143 double r0[6]={3,3,3,3,3,3};
144 double r1[6]={5,5,5,5,5,5};
145 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
146
147 if (value<0) continue;
148 if(value>maxprob) {
149 maxprob=value;
150 cout << "ir= " << ir << endl;
151// cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
152// cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
153// cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
154 cout << "MAX====> " << maxprob << endl;
155 }
156 }
157 printf("MAXprob = %.10f\n",maxprob);
158*/
159
160 p->initializePhaseSpace(getNDaug(),getDaugs());
161 EvtVector4R D1 = p->getDaug(0)->getP4();
162 EvtVector4R D2 = p->getDaug(1)->getP4();
163 EvtVector4R D3 = p->getDaug(2)->getP4();
164
165 double P1[4], P2[4], P3[4];
166 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
167 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
168 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
169
170 double value;
171 int g0[6]={2,3,1,1,2,2};
172 int spin[6]={0,0,0,2,1,1};
173 int nstates=6;
174 double r0[6]={3,3,3,3,3,3};
175 double r1[6]={5,5,5,5,5,5};
176
177 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
178 setProb(value);
179 return ;
180}
181
182//------------base---------------------------------
183void EvtDTopipi0pi0::Com_Multi(double a1[2], double a2[2], double res[2])
184{
185 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
186 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
187}
188void EvtDTopipi0pi0::Com_Divide(double a1[2], double a2[2], double res[2])
189{
190 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
191 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
192 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
193}
194double EvtDTopipi0pi0::SCADot(double a1[4], double a2[4])
195{
196 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
197 return _cal;
198}
199double EvtDTopipi0pi0::barrier(int l, double sa, double sb, double sc, double r, double mass)
200{
201 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
202 if(q < 0) q = 1e-16;
203 double z;
204 z = q*r*r;
205 double sa0;
206 sa0 = mass*mass;
207 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
208 if(q0 < 0) q0 = 1e-16;
209 double z0 = q0*r*r;
210 double F = 0.0;
211 if(l == 0) F = 1;
212 if(l == 1) F = sqrt((1+z0)/(1+z));
213 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
214 return F;
215}
216//------------------spin-------------------------------------------
217void EvtDTopipi0pi0::calt1(double daug1[4], double daug2[4], double t1[4])
218{
219 double p, pq, tmp;
220 double pa[4], qa[4];
221 for(int i=0; i<4; i++) {
222 pa[i] = daug1[i] + daug2[i];
223 qa[i] = daug1[i] - daug2[i];
224 }
225 p = SCADot(pa,pa);
226 pq = SCADot(pa,qa);
227 tmp = pq/p;
228 for(int i=0; i<4; i++) {
229 t1[i] = qa[i] - tmp*pa[i];
230 }
231}
232void EvtDTopipi0pi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
233{
234 double p, r;
235 double pa[4], t1[4];
236 calt1(daug1,daug2,t1);
237 r = SCADot(t1,t1)/3.0;
238 for(int i=0; i<4; i++) {
239 pa[i] = daug1[i] + daug2[i];
240 }
241 p = SCADot(pa,pa);
242 for(int i=0; i<4; i++) {
243 for(int j=0; j<4; j++) {
244 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
245 }
246 }
247}
248//-------------------prop--------------------------------------------
249void EvtDTopipi0pi0::propagator(double mass2, double mass, double width, double sx, double prop[2])
250{
251 double a[2], b[2];
252 a[0] = 1;
253 a[1] = 0;
254 b[0] = mass2-sx;
255 b[1] = -mass*width;
256 Com_Divide(a,b,prop);
257}
258double EvtDTopipi0pi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
259{
260 double widm = 0.;
261 double m = sqrt(sa);
262 double tmp = sb-sc;
263 double tmp1 = sa+tmp;
264 double q = 0.25*tmp1*tmp1/sa-sb;
265 if(q<0) q = 1e-16;
266 double tmp2 = mass2+tmp;
267 double q0 = 0.25*tmp2*tmp2/mass2-sb;
268 if(q0<0) q0 = 1e-16;
269 double z = q*r2;
270 double z0 = q0*r2;
271 double t = q/q0;
272 if(l == 0) {widm = sqrt(t)*mass/m;}
273 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
274 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
275 return widm;
276}
277double EvtDTopipi0pi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
278{
279 double widm = 0.;
280 double m = sqrt(sa);
281 double tmp = sb-sc;
282 double tmp1 = sa+tmp;
283 double q = 0.25*tmp1*tmp1/sa-sb;
284 if(q<0) q = 1e-16;
285 double tmp2 = mass2+tmp;
286 double q0 = 0.25*tmp2*tmp2/mass2-sb;
287 if(q0<0) q0 = 1e-16;
288 double z = q*r2;
289 double z0 = q0*r2;
290 double F = (1+z0)/(1+z);
291 double t = q/q0;
292 widm = t*sqrt(t)*mass/m*F;
293 return widm;
294}
295void EvtDTopipi0pi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
296{
297 double a[2], b[2];
298 a[0] = 1;
299 a[1] = 0;
300 b[0] = mass2-sa;
301 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
302 Com_Divide(a,b,prop);
303}
304
305void EvtDTopipi0pi0::propagatorFlatte(double mass, double width, double sa, double sb, double sc, int r, double prop[2]){
306 double q, qKsK;
307 double rhoab[2], rhoKsK[2];
308 q = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
309 if(r == 0) qKsK = 0.25*sa - 0.49368*0.49368;
310 if(r == 1) qKsK = 0.25*(sa + 3.899750596e-03)*(sa + 3.899750596e-03)/sa - 0.497614*0.497614;
311 if(q>0){
312 rhoab[0] = 2*sqrt(q/sa);
313 rhoab[1] = 0;
314 }
315 if(q<0){
316 rhoab[0] = 0;
317 rhoab[1] = 2*sqrt(-q/sa);
318 }
319 if(qKsK>0){
320 rhoKsK[0] = 2*sqrt(qKsK/sa);
321 rhoKsK[1] = 0;
322 }
323 if(qKsK<0){
324 rhoKsK[0] = 0;
325 rhoKsK[1] = 2*sqrt(-qKsK/sa);
326 }
327 double a[2], b[2];
328 a[0] = 1;
329 a[1] = 0;
330 b[0] = mass*mass - sa + 0.341*rhoab[1] + 0.892*0.341*rhoKsK[1];
331 b[1] = - (0.341*rhoab[0] + 0.892*0.341*rhoKsK[0]);
332 Com_Divide(a,b,prop);
333}
334
335void EvtDTopipi0pi0::rhoab(double sa, double sb, double sc, double res[2]) {
336 double tmp = sa+sb-sc;
337 double q = 0.25*tmp*tmp/sa-sb;
338 if(q>=0) {
339 res[0]=2.0*sqrt(q/sa);
340 res[1]=0.0;
341 } else {
342 res[0]=0.0;
343 res[1]=2.0*sqrt(-q/sa);
344 }
345}
346void EvtDTopipi0pi0::rho4Pi(double sa, double res[2]) {
347 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
348 if(temp>=0) {
349 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
350 res[1]=0.0;
351 } else {
352 res[0]=0.0;
353 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));//?????????????????????
354 }
355}
356
357void EvtDTopipi0pi0::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
358 double f = 0.5843+1.6663*sa;
359 const double M = 0.9264;//M(f0500)
360 const double mass2 = 0.85821696; // M*M
361 const double mpi2d2 = 0.00973989245;//mass_Pion2/2 = 0.0194797849/2
362 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
363 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
364 rhoab(sa,sb,sc,rho1s);
365 rhoab(mass2,sb,sc,rho1M);
366 rho4Pi(sa,rho2s);
367 rho4Pi(mass2,rho2M);
368 Com_Divide(rho1s,rho1M,rho1);
369 Com_Divide(rho2s,rho2M,rho2);
370 double a[2], b[2];
371 a[0] = 1.0;
372 a[1] = 0.0;
373 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
374 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
375 Com_Divide(a,b,prop);
376}
377void EvtDTopipi0pi0::Flatte_rhoab(double sa, double sb, double rho[2])
378{
379 double q = 1.0-(4*sb/sa);
380
381 if(q>0) {
382 rho[0]=sqrt(q);
383 rho[1]=0;
384 }
385 else if(q<0){
386 rho[0]=0;
387 rho[1]=sqrt(-q);
388 }
389}
390
391void EvtDTopipi0pi0::propagator980(double mass, double sx, double *sb, double prop[2])
392{
393 double gpipi1 = 2.0/3.0*0.165;
394 double gpipi2 = 1.0/3.0*0.165;
395 double gKK1 = 0.5*0.69465;
396 double gKK2 = 0.5*0.69465;
397
398 double unit[2]={1.0};
399 double ci[2]={0,1};
400 double rho1[2];
401 Flatte_rhoab(sx,sb[0],rho1);
402 double rho2[2];
403 Flatte_rhoab(sx,sb[1],rho2);
404 double rho3[2];
405 Flatte_rhoab(sx,sb[2],rho3);
406 double rho4[2];
407 Flatte_rhoab(sx,sb[3],rho4);
408
409 double tmp1[2]={gpipi1,0};
410 double tmp11[2];
411 double tmp2[2]={gpipi2,0};
412 double tmp22[2];
413 double tmp3[2]={gKK1,0};
414 double tmp33[2];
415 double tmp4[2]={gKK2,0};
416 double tmp44[2];
417
418 Com_Multi(tmp1,rho1,tmp11);
419 Com_Multi(tmp2,rho2,tmp22);
420 Com_Multi(tmp3,rho3,tmp33);
421 Com_Multi(tmp4,rho4,tmp44);
422
423 double tmp5[2]={tmp11[0]+tmp22[0]+tmp33[0]+tmp44[0],tmp11[1]+tmp22[1]+tmp33[1]+tmp44[1]};
424 double tmp51[2];
425 Com_Multi(tmp5, ci,tmp51);
426 double tmp6[2]={mass*mass-sx-tmp51[0], -1.0*tmp51[1]};
427
428 Com_Divide( unit,tmp6, prop);
429}
430
431//------------GS---used by rho----------------------------
432void EvtDTopipi0pi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
433{
434 double a[2], b[2];
435 double tmp = sb-sc;
436 double tmp1 = sa+tmp;
437 double q2 = 0.25*tmp1*tmp1/sa-sb;
438 if(q2<0) q2 = 1e-16;
439
440 double tmp2 = mass2+tmp;
441 double q02 = 0.25*tmp2*tmp2/mass2-sb;
442 if(q02<0) q02 = 1e-16;
443
444 double q = sqrt(q2);
445 double q0 = sqrt(q02);
446 double m = sqrt(sa);
447 double q03 = q0*q02;
448 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
449
450 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
451 double h0 = GS1*q0/mass*tmp3;
452 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
453 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
454 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
455
456 a[0] = 1.0+d*width/mass;
457 a[1] = 0.0;
458 b[0] = mass2-sa+width*f;
459 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
460 Com_Divide(a,b,prop);
461}
462
463double EvtDTopipi0pi0::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass, double rRES0, double rRES1){
464 double pR[4], pD[4];
465 double temp_PDF, v_re;
466 temp_PDF = 0.0;
467 v_re = 0.0;
468 double B[2], s1, s2, s3, sR, sD;
469 for(int i=0; i<4; i++){
470 pR[i] = P1[i] + P2[i];
471 pD[i] = pR[i] + P3[i];
472 }
473 s1 = SCADot(P1,P1);
474 s2 = SCADot(P2,P2);
475 s3 = SCADot(P3,P3);
476 sR = SCADot(pR,pR);
477 sD = SCADot(pD,pD);
478 int G[4][4];
479 for(int i=0; i!=4; i++){
480 for(int j=0; j!=4; j++){
481 if(i==j){
482 if(i==0) G[i][j] = 1;
483 else G[i][j] = -1;
484 }
485 else G[i][j] = 0;
486 }
487 }
488 if(Ang == 0){
489 B[0] = 1;
490 B[1] = 1;
491 temp_PDF = 1;
492 }
493 if(Ang == 1){
494 B[0] = barrier(1,sR,s1,s2,rRES0,mass);
495 B[1] = barrier(1,sD,sR,s3,rRES1,1.86966);
496 double t1[4], T1[4];
497 calt1(P1,P2,t1);
498 calt1(pR,P3,T1);
499 temp_PDF = 0;
500 for(int i=0; i<4; i++){
501 temp_PDF += t1[i]*T1[i]*G[i][i];
502 }
503 }
504 if(Ang == 2){
505 B[0] = barrier(2,sR,s1,s2,rRES0,mass);
506 B[1] = barrier(2,sD,sR,s3,rRES1,1.86966);
507 double t2[4][4], T2[4][4];
508 calt2(P1,P2,t2);
509 calt2(pR,P3,T2);
510 temp_PDF = 0;
511 for(int i=0; i<4; i++){
512 for(int j=0; j<4; j++){
513 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
514 }
515 }
516 }
517 v_re = temp_PDF*B[0]*B[1];
518 return v_re;
519}
520
521void 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)
522{
523 double P12[4], P23[4], P13[4];
524// double Pic[4], Pi01[4], Pi02[4], P12[4], P23[4], P13[4];
525 double cof[2], amp_PDF[2], PDF[2];
526 double spi01, spi02, scpi;
527 double s12, s13, s23;
528 for(int i=0; i<4; i++){
529 P23[i] = Pi01[i] + Pi02[i];
530 P12[i] = Pic[i] + Pi01[i];
531 P13[i] = Pic[i] + Pi02[i];
532 }
533 scpi = SCADot(Pic,Pic);
534 spi01 = SCADot(Pi01,Pi01);
535 spi02 = SCADot(Pi02,Pi02);
536 s23 = SCADot(P23,P23);
537 s12 = SCADot(P12,P12);
538 s13 = SCADot(P13,P13); //pipi0
539 double spi012[4]={mass_Pion*mass_Pion,spi02,mass_Kaon*mass_Kaon,mk0*mk0};
540// double spi022[4]={mass_Pion*mass_Pion,spi02,mass_Kaon*mass_Kaon,mk0*mk0};
541
542// double spi012[4]={spi01,spi02,mass_Kaon*mass_Kaon,mass_Kaon*mass_Kaon};
543// double spi022[4]={spi02,spi02,mass_Kaon*mass_Kaon,mass_Kaon*mass_Kaon};
544
545 double pro[2], temp_PDF, pro1[2], temp_PDF1, pro2[2], temp_PDF2, amp_tmp[2];
546 double mass1sq;
547 amp_PDF[0] = 0;
548 amp_PDF[1] = 0;
549 PDF[0] = 0;
550 PDF[1] = 0;
551 amp_tmp[0] = 0;
552 amp_tmp[1] = 0;
553 for(int i=0; i<nstates; i++) {
554 amp_tmp[0] = 0;
555 amp_tmp[1] = 0;
556 mass1sq = mass1[i]*mass1[i];
557 cof[0] = amp[i]*cos(phase[i]);
558 cof[1] = amp[i]*sin(phase[i]);
559 temp_PDF = 0;
560 if(modetype[i] == 23){//pi0pi0
561 temp_PDF = DDalitz(Pi01, Pi02, Pic, spin[i], mass1[i], r0[i],r1[i]);
562 if(g0[i]==2) propagatorsigma500(s23,spi01,spi02,pro);
563 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s23,spi01,spi02,r0[i]*r0[i],spin[i],pro);
564 if(g0[i]==3) propagator980(mass1[i],s23,spi012,pro);
565
566 if(g0[i]==0){
567 pro[0] = 1;
568 pro[1] = 0;
569 }
570 amp_tmp[0] = temp_PDF*pro[0];
571 amp_tmp[1] = temp_PDF*pro[1];
572 }
573 if(modetype[i] == 12){//pipi0
574 temp_PDF1 = DDalitz(Pic, Pi01, Pi02, spin[i], mass1[i], r0[i], r1[i]);
575 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s12,scpi,spi01,r0[i]*r0[i],spin[i],pro1);
576 if(g0[i]==2) propagatorGS(mass1sq,mass1[i],width1[i],s12,scpi,spi01,r0[i]*r0[i],pro1);
577 if(g0[i]==0){
578 pro1[0] = 1;
579 pro1[1] = 0;
580 }
581
582 temp_PDF2 = DDalitz(Pic, Pi02, Pi01, spin[i], mass1[i], r0[i], r1[i]);
583 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,scpi,spi02,r0[i]*r0[i],spin[i],pro2);
584 if(g0[i]==2) propagatorGS(mass1sq,mass1[i],width1[i],s13,scpi,spi02,r0[i]*r0[i],pro2);
585 if(g0[i]==0){
586 pro2[0] = 1;
587 pro2[1] = 0;
588 }
589 amp_tmp[0] = temp_PDF1*pro1[0] + temp_PDF2*pro2[0];
590 amp_tmp[1] = temp_PDF1*pro1[1] + temp_PDF2*pro2[1];
591
592 }
593 Com_Multi(amp_tmp,cof,amp_PDF);
594 PDF[0] += amp_PDF[0];
595 PDF[1] += amp_PDF[1];
596
597 }
598 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
599 if(value <=0) value = 1e-20;
600 Result = value;
601}
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
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
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:34
EvtDTopipi0pi0::decay
void decay(EvtParticle *p)
Definition: EvtDTopipi0pi0.cc:125
EvtDTopipi0pi0::clone
EvtDecayBase * clone()
Definition: EvtDTopipi0pi0.cc:40
EvtDTopipi0pi0::getName
void getName(std::string &name)
Definition: EvtDTopipi0pi0.cc:36
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition: EvtDecayBase.cc:533
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition: EvtDecayBase.cc:520
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cc:297
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition: EvtDecayBase.cc:504
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: EvtDecayBase.cc:482
EvtDecayProb::setProb
void setProb(double prob)
Definition: EvtDecayProb.hh:34
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:120
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:84
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1070
EvtVector4R::get
double get(int i) const
Definition: EvtVector4R.hh:179
s12
double double double double * s12
Definition: qcdloop1.h:77
s23
double double double double double * s23
Definition: qcdloop1.h:77
b
const double b
Definition: slope.cxx:9