BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
BesEvtGen-00-04-08/src/EvtGen/EvtGenModels/EvtD0ToKpipi0pi0.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: EvtD0ToKpipi0pi0.cc
11// the necessary file: EvtD0ToKpipi0pi0.hh
12//
13// Description: D0 -> K- pi+ pi0 pi0,
14// see PHYSICAL REVIEW D 99, 092008 (2019)
15//
16// Modification history:
17// Liaoyuan Dong Jan.15, 2020 Module created
18//------------------------------------------------------------------------
19//
20#include "EvtGenBase/EvtPatches.hh"
21#include "EvtGenBase/EvtParticle.hh"
22#include "EvtGenBase/EvtGenKine.hh"
23#include "EvtGenBase/EvtPDL.hh"
24#include "EvtGenBase/EvtReport.hh"
25#include "EvtGenModels/EvtD0ToKpipi0pi0.hh"
26#include "EvtGenBase/EvtComplex.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
28#include <stdlib.h>
29
30EvtD0ToKpipi0pi0::~EvtD0ToKpipi0pi0() {}
31
32void EvtD0ToKpipi0pi0::getName(std::string& model_name){
33 model_name="D0ToKpipi0pi0";
34}
35
36EvtDecayBase* EvtD0ToKpipi0pi0::clone(){
37 return new EvtD0ToKpipi0pi0;
38}
39
40void EvtD0ToKpipi0pi0::init(){
41 checkNArg(0);
42 checkNDaug(4);
43 checkSpinParent(EvtSpinType::SCALAR);
44/*
45 checkSpinDaughter(0,EvtSpinType::SCALAR);
46 checkSpinDaughter(1,EvtSpinType::SCALAR);
47 checkSpinDaughter(3,EvtSpinType::SCALAR);
48 checkSpinDaughter(4,EvtSpinType::SCALAR);
49*/
50 std::cout << "EvtD0ToKpipi0pi0 ==> Initialization !" << std::endl;
51
52 mod[0]= 1; rho[0]= 2.02; phi[0]= -0.75;
53 mod[1]= 1; rho[1]= 1.66; phi[1]= -2.90;
54 mod[2]= 0; rho[2]= 0; phi[2]= 0;
55 mod[3]= 0; rho[3]= 0; phi[3]= 0;
56 mod[4]= 0; rho[4]= 0; phi[4]= 0;
57 mod[5]= 0; rho[5]= 0; phi[5]= 0;
58 mod[6]= 0; rho[6]= 0; phi[6]= 0;
59 mod[7]= 1; rho[7]= 1; phi[7]= 0;
60 mod[8]= 1; rho[8]= 0.842; phi[8]= -2.05;
61 mod[9]= 1; rho[9]= 0.0218; phi[9]= 1.84;
62 mod[10]=0; rho[10]= 0; phi[10]= 0;
63 mod[11]=0; rho[11]= 0; phi[11]= 0;
64 mod[12]=0; rho[12]= 0; phi[12]= 0;
65 mod[13]=1; rho[13]= 0.0336; phi[13]= -1.55;
66 mod[14]=1; rho[14]= 0.109; phi[14]= -1.35;
67 mod[15]=1; rho[15]= 0.196; phi[15]= -2.07;
68 mod[16]=0; rho[16]= 0; phi[16]= 0;
69 mod[17]=0; rho[17]= 0; phi[17]= 0;
70 mod[18]=0; rho[18]= 0; phi[18]= 0;
71 mod[19]=1; rho[19]= 0.363; phi[19]= 1.93;
72 mod[20]=0; rho[20]= 0; phi[20]= 0;
73 mod[21]=0; rho[21]= 0; phi[21]= 0;
74 mod[22]=0; rho[22]= 0; phi[22]= 0;
75 mod[23]=1; rho[23]= 0.555; phi[23]= 0.44;
76 mod[24]=1; rho[24]= 0.526; phi[24]= -1.84;
77 mod[25]=0; rho[25]= 0; phi[25]= 0;
78 mod[26]=1; rho[26]= 1; phi[26]= 0.64;
79 mod[27]=0; rho[27]= 0; phi[27]= 0;
80 mod[28]=0; rho[28]= 0; phi[28]= 0;
81 mod[29]=1; rho[29]= 3.34; phi[29]= -0.02;
82 mod[30]=0; rho[30]= 0; phi[30]= 0;
83 mod[31]=0; rho[31]= 0; phi[31]= 0;
84 mod[32]=0; rho[32]= 0; phi[32]= 0;
85 mod[33]=0; rho[33]= 0; phi[33]= 0;
86 mod[34]=1; rho[34]= 1.76; phi[34]= -2.39;
87 mod[35]=1; rho[35]= 0.175; phi[35]= 1.59;
88 mod[36]=1; rho[36]= 0.397; phi[36]= 1.45;
89 mod[37]=0; rho[37]= 0; phi[37]= 0;
90 mod[38]=0; rho[38]= 0; phi[38]= 0;
91 mod[39]=0; rho[39]= 0; phi[39]= 0;
92 mod[40]=0; rho[40]= 0; phi[40]= 0;
93 mod[41]=1; rho[41]= 1.02; phi[41]= 0.52;
94 mod[42]=0; rho[42]= 0; phi[42]= 0;
95 mod[43]=0; rho[43]= 0; phi[43]= 0;
96 mod[44]=0; rho[44]= 0; phi[44]= 0;
97 mod[45]=0; rho[45]= 0; phi[45]= 0;
98 mod[46]=1; rho[46]= 0.146; phi[46]= 1.24;
99 mod[47]=1; rho[47]= 0.0978; phi[47]= -2.89;
100 mod[48]=1; rho[48]= 0.233; phi[48]= 2.41;
101 mod[49]=0; rho[49]= 0; phi[49]= 0;
102 mod[50]=1; rho[50]= 0.424; phi[50]= -0.94;
103 mod[51]=1; rho[51]= 1.03; phi[51]= -1.93;
104 mod[52]=0; rho[52]= 0; phi[52]= 0;
105 mod[53]=0; rho[53]= 0; phi[53]= 0;
106 mod[54]=1; rho[54]= 0.474; phi[54]= -1.17;
107 mod[55]=0; rho[55]= 0; phi[55]= 0;
108 mod[56]=0; rho[56]= 0; phi[56]= 0;
109 mod[57]=0; rho[57]= 0; phi[57]= 0;
110 mod[58]=0; rho[58]= 0; phi[58]= 0;
111 mod[59]=0; rho[59]= 0; phi[59]= 0;
112 mod[60]=0; rho[60]= 0; phi[60]= 0;
113 mod[61]=1; rho[61]= 6.74; phi[61]= -1.74;
114 mod[62]=0; rho[62]= 0; phi[62]= 0;
115 mod[63]=0; rho[63]= 0; phi[63]= 0;
116 mod[64]=0; rho[64]= 0; phi[64]= 0;
117 mod[65]=0; rho[65]= 0; phi[65]= 0;
118 mod[66]=1; rho[66]= 1.54; phi[66]= -2.93;
119 mod[67]=1; rho[67]= 1.36; phi[67]= 2.23;
120
121 mass[0]= 1.23; width[0]= 0.50204;
122 mass[1]= 1.2723; width[1]= 0.09;
123 mass[2]= 0.89166; width[2]= 0.0508;
124 mass[3]= 0.89581; width[3]= 0.0474;
125 mass[4]= 0.77511; width[4]= 0.1491;
126
127 for (int i=0; i<69; i++) {
128 cout << i << "rho,phi = " << rho[i] << ", "<< phi[i] << endl;
129 }
130 mD = 1.86486;
131 rRes = 3.0;
132 rD = 5.0;
133 metap = 0.95778;
134 mk0 = 0.497614;
135 mass_Kaon = 0.49368;
136 mass_Pion = 0.13957;
137 math_pi = 3.1415926;
138 mass_Pi0 = 0.1349766;
139 mkstrm = 0.89594;
140 mkstr0 = 0.89594;
141
142 pi = 3.1415926;
143 mpi = 0.13957;
144 g1 = 0.5468;
145 g2 = 0.23;
146
147 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
148 int EE[4][4][4][4] =
149 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
150 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
151 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
152 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
153 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
154 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
155 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
156 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
157 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
158 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
159 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
160 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
161 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
162 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
163 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
164 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
165 for (int i=0; i<4; i++) {
166 for (int j=0; j<4; j++) {
167 G[i][j] = GG[i][j];
168 for (int k=0; k<4; k++) {
169 for (int l=0; l<4; l++) {
170 E[i][j][k][l] = EE[i][j][k][l];
171 }
172 }
173 }
174 }
175}
176
177void EvtD0ToKpipi0pi0::initProbMax() {
178 setProbMax(8060.0);
179}
180
181void EvtD0ToKpipi0pi0::decay( EvtParticle *p){
182/*
183 double maxprob = 0.0;
184 for(int ir=0;ir<=60000000;ir++){
185 p->initializePhaseSpace(getNDaug(),getDaugs());
186 EvtVector4R Km0 = p->getDaug(0)->getP4();
187 EvtVector4R pi1 = p->getDaug(1)->getP4();
188 EvtVector4R pi2 = p->getDaug(2)->getP4();
189 EvtVector4R pi3 = p->getDaug(3)->getP4();
190 double Km[4],Pip[4],Pi01[4],Pi02[4];
191 Km[0] = Km0.get(0); Pip[0] = pi1.get(0); Pi01[0] = pi2.get(0); Pi02[0] = pi3.get(0);
192 Km[1] = Km0.get(1); Pip[1] = pi1.get(1); Pi01[1] = pi2.get(1); Pi02[1] = pi3.get(1);
193 Km[2] = Km0.get(2); Pip[2] = pi1.get(2); Pi01[2] = pi2.get(2); Pi02[2] = pi3.get(2);
194 Km[3] = Km0.get(3); Pip[3] = pi1.get(3); Pi01[3] = pi2.get(3); Pi02[3] = pi3.get(3);
195 double Prob = calPDF(Km, Pip, Pi01, Pi02);
196 if(Prob>maxprob) {
197 maxprob=Prob;
198 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
199 }
200 }
201 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
202*/
203 p->initializePhaseSpace(getNDaug(),getDaugs());
204 EvtVector4R Km0 = p->getDaug(0)->getP4();
205 EvtVector4R pi1 = p->getDaug(1)->getP4();
206 EvtVector4R pi2 = p->getDaug(2)->getP4();
207 EvtVector4R pi3 = p->getDaug(3)->getP4();
208
209 double Km[4],Pip[4],Pi01[4],Pi02[4];
210 Km[0] = Km0.get(0); Pip[0] = pi1.get(0); Pi01[0] = pi2.get(0); Pi02[0] = pi3.get(0);
211 Km[1] = Km0.get(1); Pip[1] = pi1.get(1); Pi01[1] = pi2.get(1); Pi02[1] = pi3.get(1);
212 Km[2] = Km0.get(2); Pip[2] = pi1.get(2); Pi01[2] = pi2.get(2); Pi02[2] = pi3.get(2);
213 Km[3] = Km0.get(3); Pip[3] = pi1.get(3); Pi01[3] = pi2.get(3); Pi02[3] = pi3.get(3);
214 double prob = calPDF(Km, Pip, Pi01, Pi02);
215 setProb(prob);
216 return;
217}
218
219double EvtD0ToKpipi0pi0::calPDF(double Km[], double Pip[], double Pi01[], double Pi02[])
220{
221 Km[0] = sqrt(mass_Kaon*mass_Kaon + Km[1]*Km[1] + Km[2]*Km[2] + Km[3]*Km[3]);
222 Pip[0] = sqrt(mass_Pion*mass_Pion + Pip[1]*Pip[1] + Pip[2]*Pip[2] + Pip[3]*Pip[3]);
223 Pi01[0] = sqrt(mass_Pi0*mass_Pi0 + Pi01[1]*Pi01[1] + Pi01[2]*Pi01[2] + Pi01[3]*Pi01[3]);
224 Pi02[0] = sqrt(mass_Pi0*mass_Pi0 + Pi02[1]*Pi02[1] + Pi02[2]*Pi02[2] + Pi02[3]*Pi02[3]);
225
226 EvtComplex PDF[68];
227 int g[3];
228 //-------PHSP---------
229 PDF[0] = PHSP(Km,Pip) + PHSP(Km,Pip);
230 PDF[1] = PHSP(Km,Pi01)+PHSP(Km,Pi02);
231 //-------D2PP,P2VP------
232// PDF[2] = D2PP_P2VP(Pi01,Pi02,Km,Pip,0) + D2PP_P2VP(Pi02,Pi01,Km,Pip,0);
233// PDF[3] = D2PP_P2VP(Pi01,Pip,Km,Pi02,10) + D2PP_P2VP(Pi02,Pip,Km,Pi01,10);
234// PDF[4] = D2PP_P2VP(Pip,Pi01,Km,Pi02,20) + D2PP_P2VP(Pip,Pi02,Km,Pi01,20);
235// PDF[5] = D2PP_P2VP(Pi01,Km,Pip,Pi02,1) + D2PP_P2VP(Pi02,Km,Pip,Pi01,1);
236// PDF[6] = D2PP_P2VP(Km,Pi01,Pip,Pi02,11) + D2PP_P2VP(Km,Pi02,Pip,Pi01,11);
237 //----------D2AP,A2VP--------------
238 g[0] = 1; g[1] = 1; g[2] = 0;
239 PDF[7] = D2AP_A2VP(Km,Pi01,Pip,Pi02,g,0) + D2AP_A2VP(Km,Pi02,Pip,Pi01,g,0);
240 g[2] = 2;
241 PDF[8] = D2AP_A2VP(Km,Pi01,Pip,Pi02,g,0) + D2AP_A2VP(Km,Pi02,Pip,Pi01,g,0);
242 g[2] = 0;
243 PDF[9] = D2AP_A2VP(Pip,Pi01,Km,Pi02,g,1) + D2AP_A2VP(Pip,Pi02,Km,Pi01,g,1);
244// g[2] = 2;
245// PDF[10] = D2AP_A2VP(Pip,Pi01,Km,Pi02,g,1) + D2AP_A2VP(Pip,Pi02,Km,Pi01,g,1);
246// g[2] = 0;
247// PDF[11] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
248// g[2] = 2;
249// PDF[12] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
250 g[2] = 0;
251 PDF[13] = D2AP_A2VP(Pi01,Pi02,Km,Pip,g,31) + D2AP_A2VP(Pi02,Pi01,Km,Pip,g,31);
252 g[2] = 2;
253 PDF[14] = D2AP_A2VP(Pi01,Pi02,Km,Pip,g,31) + D2AP_A2VP(Pi02,Pi01,Km,Pip,g,31);
254 g[2] = 0;
255 PDF[15] = D2AP_A2VP(Pi01,Km,Pip,Pi02,g,3) + D2AP_A2VP(Pi02,Km,Pip,Pi01,g,3);
256// g[2] = 2;
257// PDF[16] = D2AP_A2VP(Pi01,Km,Pip,Pi02,g,3) + D2AP_A2VP(Pi02,Km,Pip,Pi01,g,3);
258 g[0] = 1; g[1] = 0; g[2] = 0;
259// PDF[17] = D2AP_A2VP(Km,Pi01,Pip,Pi02,g,0) + D2AP_A2VP(Km,Pi02,Pip,Pi01,g,0);
260// g[2] = 2;
261// PDF[18] = D2AP_A2VP(Km,Pi01,Pip,Pi02,g,0) + D2AP_A2VP(Km,Pi02,Pip,Pi01,g,0);
262 g[2] = 0;
263 PDF[19] = D2AP_A2VP(Pip,Pi01,Km,Pi02,g,1) + D2AP_A2VP(Pip,Pi02,Km,Pi01,g,1);
264// g[2] = 2;
265// PDF[20] = D2AP_A2VP(Pip,Pi01,Km,Pi02,g,1) + D2AP_A2VP(Pip,Pi02,Km,Pi01,g,1);
266// g[2] = 0;
267// PDF[21] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
268// g[2] = 2;
269// PDF[22] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
270 g[2] = 0;
271 PDF[23] = D2AP_A2VP(Pi01,Pi02,Km,Pip,g,31) + D2AP_A2VP(Pi02,Pi01,Km,Pip,g,31);
272 g[2] = 2;
273 PDF[24] = D2AP_A2VP(Pi01,Pi02,Km,Pip,g,31) + D2AP_A2VP(Pi02,Pi01,Km,Pip,g,31);
274// g[2] = 0;
275// PDF[25] = D2AP_A2VP(Pi01,Km,Pip,Pi02,g,3) + D2AP_A2VP(Pi02,Km,Pip,Pi01,g,3);
276 g[2] = 2;
277 PDF[26] = D2AP_A2VP(Pi01,Km,Pip,Pi02,g,3) + D2AP_A2VP(Pi02,Km,Pip,Pi01,g,3);
278 //--------D2AP,A2SP-----------------------------------
279// PDF[27] = D2AP_A2SP(Km,Pi01,Pip,Pi02,0) + D2AP_A2SP(Km,Pi02,Pip,Pi01,0);
280// PDF[28] = D2AP_A2SP(Km,Pip,Pi01,Pi02,10) + D2AP_A2SP(Km,Pip,Pi02,Pi01,10);
281 PDF[29] = D2AP_A2SP(Pi01,Pi02,Km,Pip,1) + D2AP_A2SP(Pi02,Pi01,Km,Pip,1);
282// PDF[30] = D2AP_A2SP(Pi01,Pip,Km,Pi02,11) + D2AP_A2SP(Pi02,Pip,Km,Pi01,11);
283// PDF[31] = D2AP_A2SP(Pip,Pi01,Km,Pi02,21) + D2AP_A2SP(Pip,Pi02,Km,Pi01,21);
284// PDF[32] = D2AP_A2SP(Pi01,Km,Pip,Pi02,31) + D2AP_A2SP(Pi02,Km,Pip,Pi01,31);
285// PDF[33] = D2AP_A2SP(Pip,Km,Pi01,Pi02,41) + D2AP_A2SP(Pip,Km,Pi02,Pi01,41);
286 //--------D2VS----------------------------
287 PDF[34] = D2VS(Pip,Pi01,Km,Pi02,1,0) + D2VS(Pip,Pi02,Km,Pi01,1,0);
288 PDF[35] = D2VS(Km,Pi01,Pip,Pi02,1,1) + D2VS(Km,Pi02,Pip,Pi01,1,1);
289 PDF[36] = D2VS(Km,Pip,Pi01,Pi02,1,11) + D2VS(Km,Pip,Pi02,Pi01,1,11);
290// PDF[37] = D2VS(Pi01,Pi02,Km,Pip,0,10) + D2VS(Pi02,Pi01,Km,Pip,0,10);
291// PDF[38] = D2VS(Pip,Pi01,Km,Pi02,0,0) + D2VS(Pip,Pi02,Km,Pi01,0,0);
292// PDF[39] = D2VS(Km,Pip,Pi01,Pi02,0,11) + D2VS(Km,Pip,Pi02,Pi01,0,11);
293// PDF[40] = D2VS(Km,Pi01,Pip,Pi02,0,1) + D2VS(Km,Pi02,Pip,Pi01,0,1);
294 //---------D2VP,V2VP----------------------
295 PDF[41] = D2VP_V2VP(Pi01,Pip,Km,Pi02,0) + D2VP_V2VP(Pi02,Pip,Km,Pi01,0);
296// PDF[42] = D2VP_V2VP(Pip,Pi01,Km,Pi02,10) + D2VP_V2VP(Pip,Pi02,Km,Pi01,10);
297// PDF[43] = D2VP_V2VP(Pi01,Pi02,Km,Pip,20) + D2VP_V2VP(Pi02,Pi01,Km,Pip,20);
298// PDF[44] = D2VP_V2VP(Pi01,Km,Pip,Pi02,1) + D2VP_V2VP(Pi02,Km,Pip,Pi01,1);
299// PDF[45] = D2VP_V2VP(Km,Pi01,Pip,Pi02,2) + D2VP_V2VP(Km,Pi02,Pip,Pi01,2);
300 //-----------D2VV--------------------------
301 g[0] = 1; g[1] = 1; g[2] = 0;
302 PDF[46] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
303 g[2] = 1;
304 PDF[47] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
305 g[2] = 2;
306 PDF[48] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
307 g[0] = 0; g[1] = 1; g[2] = 0;
308// PDF[49] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
309 g[2] = 1;
310 PDF[50] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
311 g[2] = 2;
312 PDF[51] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
313 g[0] = 1; g[1] = 0; g[2] = 0;
314// PDF[52] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
315// g[2] = 1;
316// PDF[53] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
317 g[2] = 2;
318 PDF[54] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
319 g[0] = 1; g[1] = 0; g[2] = 0;
320// PDF[55] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
321// g[2] = 1;
322// PDF[56] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
323// g[2] = 2;
324// PDF[57] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
325 g[0] = 0; g[1] = 0; g[2] = 0;
326// PDF[58] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
327// g[2] = 1;
328// PDF[59] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
329// g[2] = 2;
330// PDF[60] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
331 g[0] = 0; g[1] = 0; g[2] = 0;
332 PDF[61] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
333 g[2] = 1;
334// PDF[62] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
335 g[2] = 2;
336// PDF[63] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
337 //----------D2TS--------------------
338// PDF[64] = D2TS(Km,Pip,Pi01,Pi02,0) + D2TS(Km,Pip,Pi02,Pi01,0);
339// PDF[65] = D2TS(Km,Pi01,Pip,Pi02,10) + D2TS(Km,Pi02,Pip,Pi01,10);
340 PDF[66] = D2TS(Pi02,Pi01,Km,Pip,1) + D2TS(Pi01,Pi02,Km,Pip,1);
341 PDF[67] = D2TS(Pip,Pi01,Km,Pi02,11) + D2TS(Pip,Pi02,Km,Pi01,11);
342
343//------------------------------------------
344 EvtComplex cof;
345 EvtComplex pdf, module;
346 pdf = EvtComplex(0,0);
347 for(int i=0; i < 68; i++){
348 if (mod[i]==0) continue;
349 cof = EvtComplex(rho[i]*cos(phi[i]),rho[i]*sin(phi[i]));
350 pdf = pdf + cof*PDF[i];
351 }
352 module = conj(pdf)*pdf;
353 double value;
354 value = real(module);
355 return (value <= 0) ? 1e-20 : value;
356}
357EvtComplex EvtD0ToKpipi0pi0::KPiSFormfactor(double sa, double sb, double sc, double r)
358{
359 double m1430 = 1.463;
360 double sa0 = m1430*m1430;
361 double w1430 = 0.233;
362 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
363 if(q0<0) q0 = 1e-16;
364 double qs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
365 double q = sqrt(qs);
366 double width = w1430*q*m1430/sqrt(sa*q0);
367 double temp_R = atan(m1430*width/(sa0-sa));
368 if(temp_R<0) temp_R += math_pi;
369 double deltaR = -5.31 + temp_R;
370 double temp_F = atan(2*1.07*q/(2+(-1.8)*1.07*qs));
371 if(temp_F<0) temp_F += math_pi;
372 double deltaF = 2.33 + temp_F;
373 EvtComplex cR(cos(deltaR), sin(deltaR));
374 EvtComplex cF(cos(deltaF), sin(deltaF));
375 EvtComplex amp = 0.8*sin(deltaF)*cF + sin(deltaR)*cR*cF*cF;
376 return amp;
377}
378EvtComplex EvtD0ToKpipi0pi0::D2VV(double P1[], double P2[], double P3[], double P4[], int g[], int flag)
379{
380 double t1V1[4], t1V2[4], t1D[4], t2D[4][4];
381 double temp_PDF = 0;
382 EvtComplex amp_PDF(0,0);
383 EvtComplex pro[2];
384
385 double sa[3], sb[3], sc[3], B[3];
386 double pV1[4], pV2[4], pD[4];
387 for(int i=0; i!=4; i++){
388 pV1[i] = P1[i] + P2[i];
389 pV2[i] = P3[i] + P4[i];
390 pD[i] = pV1[i] + pV2[i];
391 }
392 sa[0] = dot(pV1,pV1);
393 sb[0] = dot(P1,P1);
394 sc[0] = dot(P2,P2);
395 sa[1] = dot(pV2,pV2);
396 sb[1] = dot(P3,P3);
397 sc[1] = dot(P4,P4);
398 sa[2] = dot(pD,pD);
399 sb[2] = sa[0];
400 sc[2] = sa[1];
401 if(g[0] == 1){
402 if(flag == 0)pro[0] = propagatorRBW(mass[2],width[2],sa[0],sb[0],sc[0],rRes,1);
403 if(flag == 1)pro[0] = propagatorRBW(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
404 }
405 if(g[1] == 1){
406 if(flag == 0) pro[1] = propagatorGS(mass[4],width[4],sa[1],sb[1],sc[1],rRes,1);//rho
407 if(flag == 1) pro[1] = 1;
408 }
409 if(g[0] == 0) pro[0] = 1;
410 if(g[1] == 0) pro[1] = 1;
411 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
412 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
413 calt1(P1,P2,t1V1);
414 calt1(P3,P4,t1V2);
415 if(g[2] == 0){
416 for(int i=0; i!=4; i++){
417 temp_PDF += (G[i][i])*t1V1[i]*t1V2[i];
418 }
419 B[2] = 1;
420 }
421 if(g[2] == 1){
422 calt1(pV1,pV2,t1D);
423 for(int i=0; i!=4; i++){
424 for(int j=0; j!=4; j++){
425 for(int k=0; k!=4; k++){
426 for(int l=0; l!=4; l++){
427 temp_PDF += E[i][j][k][l]*pD[i]*t1D[j]*t1V1[k]*t1V2[l]*(G[i][i])*(G[j][j])*(G[l][l])*(G[k][k]);
428 }
429 }
430 }
431 }
432 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
433 }
434 if(g[2] == 2){
435 calt2(pV1,pV2,t2D);
436 for(int i=0; i!=4; i++){
437 for(int j=0; j!=4; j++){
438 temp_PDF += t2D[i][j]*t1V1[i]*t1V2[j]*(G[i][i])*(G[j][j]);
439 }
440 }
441 B[2] = barrier(2,sa[2],sb[2],sc[2],rD);
442 }
443 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro[0]*pro[1];
444 return amp_PDF;
445}
446
447EvtComplex EvtD0ToKpipi0pi0::D2AP_A2VP(double P1[], double P2[], double P3[], double P4[], int g[], int flag)
448{
449 double temp_PDF = 0;
450 EvtComplex amp_PDF(0,0);
451 EvtComplex pro[2];
452 double t1V[4], t1D[4], t2A[4][4];
453 double sa[3], sb[3], sc[3], B[3];
454 double pV[4],pA[4],pD[4];
455 for(int i=0; i!=4; i++){
456 pV[i] = P3[i] + P4[i];
457 pA[i] = pV[i] + P2[i];
458 pD[i] = pA[i] + P1[i];
459 }
460 sa[0] = dot(pV,pV);
461 sb[0] = dot(P3,P3);
462 sc[0] = dot(P4,P4);
463 sa[1] = dot(pA,pA);
464 sb[1] = sa[0];
465 sc[1] = dot(P2,P2);
466 sa[2] = dot(pD,pD);
467 sb[2] = sa[1];
468 sc[2] = dot(P1,P1);
469 if(g[0] == 1){
470 if(flag == 0 || flag == 3) pro[0] = propagatorGS(mass[4],width[4],sa[0],sb[0],sc[0],rRes,1);
471 else if(flag == 1 || flag == 21) pro[0] = propagatorRBW(mass[2],width[2],sa[0],sb[0],sc[0],rRes,1);
472 else if(flag == 31) pro[0] = propagatorRBW(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
473 }
474 else if(g[0] == 0) pro[0] = 1;
475 if(g[1] == 1){
476 if(flag == 0) pro[1] = propagatorRBW(mass[0],width[0],sa[1],sb[1],sc[1],rRes,g[2]);
477 if(flag == 1 || flag == 21 || flag == 31 || flag == 3) pro[1] = propagatorRBW(mass[1],width[1],sa[1],sb[1],sc[1],rRes,g[2]);
478 }
479 else if(g[1] == 0) pro[1] = 1;
480 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
481 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
482 calt1(P3,P4,t1V);
483 calt1(pA,P1,t1D);
484 if(g[2] == 0){
485 for(int i=0; i!=4; i++){
486 for(int j=0; j!=4; j++){
487 temp_PDF += t1D[i]*(G[i][j] - pA[i]*pA[j]/sa[1])*t1V[j]*(G[i][i])*(G[j][j]);
488 }
489 }
490 B[1] = 1;
491 }
492 else if(g[2] == 2){
493 calt2(pV,P2,t2A);
494 for(int i=0; i!=4; i++){
495 for(int j=0; j!=4; j++){
496 temp_PDF += t1D[i]*t2A[i][j]*t1V[j]*(G[i][i])*(G[j][j]);
497 }
498 }
499 B[1] = barrier(2,sa[1],sb[1],sc[1],rRes);
500 }
501 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro[0]*pro[1];
502 return amp_PDF;
503}
504
505EvtComplex EvtD0ToKpipi0pi0::D2AP_A2SP(double P1[], double P2[], double P3[], double P4[], int flag)
506{
507 //flag = 0, S = rho; flag = 1, S = K*
508 double temp_PDF = 0;
509 EvtComplex amp_PDF(0,0);
510 EvtComplex pro;
511 double sa[3], sb[3], sc[3], B[3];
512 double t1D[4], t1A[4];
513 double pS[4], pA[4], pD[4];
514 for(int i=0; i!=4; i++){
515 pS[i] = P3[i] + P4[i];
516 pA[i] = pS[i] + P2[i];
517 pD[i] = pA[i] + P1[i];
518 }
519 sa[0] = dot(pS,pS);
520 sb[0] = dot(P3,P3);
521 sc[0] = dot(P4,P4);
522 sa[1] = dot(pA,pA);
523 sb[1] = sa[0];
524 sc[1] = dot(P2,P2);
525 sa[2] = dot(pD,pD);
526 sb[2] = sa[1];
527 sc[2] = dot(P1,P1);
528 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
529 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
530 calt1(pA,P1,t1D);
531 calt1(pS,P2,t1A);
532 for(int i=0; i!=4; i++){
533 temp_PDF += t1D[i]*t1A[i]*(G[i][i]);
534 }
535 amp_PDF = temp_PDF*B[1]*B[2];
536 if(flag == 1 || flag == 11 || flag == 21 ) amp_PDF = amp_PDF*KPiSFormfactor(sa[0],sb[0],sc[0],rRes);
537 return amp_PDF;
538}
539
540EvtComplex EvtD0ToKpipi0pi0::D2PP_P2VP(double P1[], double P2[], double P3[], double P4[], int flag)
541{
542 //modeidx = 0 :(K*0 pi0)pi0
543 //modeidx = 10:(K*- pi+)pi0
544 //modeidx = 20:(K*- pi0)pi+
545 //modeidx = 1 :(K- rho+)pi0
546 //modeidx = 11:K-(rho+ pi0)
547 double temp_PDF = 0;
548 EvtComplex amp(0,0);
549 EvtComplex prop;
550 double sa[3], sb[3], sc[3], B[3];
551 double t1V[4];
552 double pV[4], pP[4], pD[4];
553 for(int i=0; i!=4; i++){
554 pV[i] = P3[i] + P4[i];
555 pP[i] = pV[i] + P2[i];
556 pD[i] = pP[i] + P1[i];
557 }
558 sa[0] = dot(pV,pV);
559 sb[0] = dot(P3,P3);
560 sc[0] = dot(P4,P4);
561 sa[1] = dot(pP,pP);
562 sb[1] = sa[0];
563 sc[1] = dot(P2,P2);
564 sa[2] = dot(pD,pD);
565 sb[2] = sa[1];
566 sc[2] = dot(P1,P1);
567 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
568 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
569 if(flag == 0) prop = propagatorRBW(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
570 else if(flag == 10 || 20) prop = propagatorRBW(mass[2],width[2],sa[0],sb[0],sc[0],rRes,1);
571 else if(flag == 1 || 11) prop = propagatorGS(mass[4],width[4],sa[0],sb[0],sc[0],rRes,1);
572 calt1(P3,P4,t1V);
573 for(int i=0; i!=4; i++){
574 temp_PDF += P2[i]*t1V[i]*(G[i][i]);
575 }
576 amp = temp_PDF*B[0]*B[1]*prop;
577 return amp;
578}
579
580EvtComplex EvtD0ToKpipi0pi0::D2VP_V2VP(double P1[], double P2[], double P3[], double P4[], int flag)
581{
582 double temp_PDF = 0;
583 EvtComplex amp_PDF(0,0);
584 EvtComplex pro;
585 double sa[3], sb[3], sc[3], B[3];
586 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
587 for(int i=0; i!=4; i++){
588 pV2[i] = P3[i] + P4[i];
589 qV2[i] = P3[i] - P4[i];
590 pV1[i] = pV2[i] + P2[i];
591 qV1[i] = pV2[i] - P2[i];
592 pD[i] = pV1[i] + P1[i];
593 }
594 for(int i=0; i!=4; i++){
595 for(int j=0; j!=4; j++){
596 for(int k=0; k!=4; k++){
597 for(int l=0; l!=4; l++){
598 temp_PDF += E[i][j][k][l]*pV1[i]*qV1[j]*P1[k]*qV2[l]*(G[i][i])*(G[j][j])*(G[k][k])*(G[l][l]);
599 }
600 }
601 }
602 }
603 sa[0] = dot(pV2,pV2);
604 sb[0] = dot(P3,P3);
605 sc[0] = dot(P4,P4);
606 sa[1] = dot(pV1,pV1);
607 sb[1] = sa[0];
608 sc[1] = dot(P2,P2);
609 sa[2] = dot(pD,pD);
610 sb[2] = sa[1];
611 sc[2] = dot(P1,P1);
612 if(flag == 0 || flag == 10) pro = propagatorRBW(mass[2],width[2],sa[0],sb[0],sc[0],rRes,1);
613 else if(flag == 20) pro = propagatorRBW(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
614 else if(flag == 1 || flag == 2) pro = propagatorGS(mass[4],width[4],sa[0],sb[0],sc[0],rRes,1);
615 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
616 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
617 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
618 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro;
619 return amp_PDF;
620}
621
622EvtComplex EvtD0ToKpipi0pi0::D2VS(double P1[], double P2[], double P3[], double P4[], int g, int flag)
623{
624 double temp_PDF = 0;
625 EvtComplex amp_PDF(0,0);
626 EvtComplex pro;
627 double sa[3], sb[3], sc[3], B[3];
628 double t1D[4], t1V[4];
629 double pS[4], pV[4], pD[4];
630 for(int i=0; i!=4; i++){
631 pS[i] = P3[i] + P4[i];
632 pV[i] = P1[i] + P2[i];
633 pD[i] = pS[i] + pV[i];
634 }
635 sa[0] = dot(pS,pS);
636 sb[0] = dot(P3,P3);
637 sc[0] = dot(P4,P4);
638 sa[1] = dot(pV,pV);
639 sb[1] = dot(P1,P1);
640 sc[1] = dot(P2,P2);
641 sa[2] = dot(pD,pD);
642 sb[2] = sa[0];
643 sc[2] = sa[1];
644 if(g == 1){
645 if(flag == 0 ) pro = propagatorGS(mass[4],width[4],sa[1],sb[1],sc[1],rRes,1);
646 else if(flag == 1) pro = propagatorRBW(mass[2],width[2],sa[1],sb[1],sc[1],rRes,1);
647 else if(flag == 11) pro = propagatorRBW(mass[3],width[3],sa[1],sb[1],sc[1],rRes,1);
648 else if(flag == 10 ) pro = 1;
649 }
650 else if(g == 0) pro = 1;
651 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
652 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
653 calt1(P1,P2,t1V);
654 calt1(pS,pV,t1D);
655 for(int i=0; i!=4; i++){
656 temp_PDF += G[i][i]*t1D[i]*t1V[i];
657 }
658 amp_PDF = temp_PDF*B[1]*B[2]*pro;
659 if(flag == 0 || flag == 10) amp_PDF *= KPiSFormfactor(sa[0],sb[0],sc[0],rRes);
660 // if(modeidx == 1 || modeidx == 11) amp_PDF *= -1.0; /// why ???????
661 return amp_PDF;
662}
663
664EvtComplex EvtD0ToKpipi0pi0::D2TS(double P1[], double P2[], double P3[], double P4[], int flag)
665{
666 // flag == 0 KPiT. 1 PiPiT
667 double temp_PDF = 0;
668 EvtComplex amp_PDF(0,0);
669 double sa[3], sb[3], sc[3], B[3];
670 double t2D[4][4], t2T[4][4];
671 double pS[4], pT[4], pD[4];
672 for(int i=0; i!=4; i++){
673 pS[i] = P3[i] + P4[i];
674 pT[i] = P1[i] + P2[i];
675 pD[i] = pT[i] + pS[i];
676 }
677 sa[0] = dot(pT,pT);
678 sb[0] = dot(P1,P1);
679 sc[0] = dot(P2,P2);
680 sa[1] = dot(pS,pS);
681 sb[1] = dot(P3,P3);
682 sc[1] = dot(P4,P4);
683 sa[2] = dot(pD,pD);
684 sb[2] = sa[0];
685 sc[2] = sa[1];
686 B[0] = barrier(2,sa[0],sb[0],sc[0],rRes);
687 B[2] = barrier(2,sa[2],sb[2],sc[2],rD);
688 calt2(P1,P2,t2T);
689 calt2(pT,pS,t2D);
690 for(int i=0; i!=4; i++){
691 for(int j=0; j!=4; j++){
692 temp_PDF += t2D[i][j]*t2T[j][i]*(G[i][i])*(G[j][j]);
693 }
694 }
695 amp_PDF = temp_PDF*B[0]*B[2];
696 if(flag == 1 || flag == 11){ amp_PDF = amp_PDF*KPiSFormfactor(sa[1],sb[1],sc[1],rRes);}
697 return amp_PDF;
698}
699
700EvtComplex EvtD0ToKpipi0pi0::PHSP(double P1[], double P2[])
701{
702 EvtComplex amp_PDF(0,0);
703 double sa,sb,sc;
704 double KPi[4];
705 for(int i=0; i!=4; i++){
706 KPi[i] = P1[i] + P2[i];
707 }
708 sa = dot(KPi,KPi);
709 sb = dot(P1,P1);
710 sc = dot(P2,P2);
711 amp_PDF = KPiSFormfactor(sa,sb,sc,rRes);
712 return amp_PDF;
713}
714
715EvtComplex EvtD0ToKpipi0pi0::propogator(double mass, double width, double sx)const
716{
717 EvtComplex ci(0,1);
718 EvtComplex prop=1.0/(mass*mass-sx-ci*mass*width);
719 return prop;
720}
721double EvtD0ToKpipi0pi0::wid(double mass, double sa, double sb, double sc, double r, int l)const
722{
723 double widm(0.), q(0.), q0(0.);
724 double sa0 = mass*mass;
725 double m = sqrt(sa);
726 q = Qabcs(sa,sb,sc);
727 q0 = Qabcs(sa0,sb,sc);
728 double z,z0;
729 z = q*r*r;
730 z0 = q0*r*r;
731 double t = q/q0;
732 double F(0.);
733 if(l == 0) F = 1;
734 if(l == 1) F = sqrt((1+z0)/(1+z));
735 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
736 widm = pow(t,l+0.5)*mass/m*F*F;
737 return widm;
738}
739double EvtD0ToKpipi0pi0::h(double m, double q)const
740{
741 double h(0.);
742 h = 2/pi*q/m*log((m+2*q)/(2*mpi));
743 return h;
744}
745double EvtD0ToKpipi0pi0::dh(double mass, double q0)const
746{
747 double dh = h(mass,q0)*(1.0/(8*q0*q0)-1.0/(2*mass*mass))+1.0/(2*pi*mass*mass);
748 return dh;
749}
750double EvtD0ToKpipi0pi0::f(double mass, double sx, double q0, double q)const
751{
752 double m = sqrt(sx);
753 double f = mass*mass/(pow(q0,3))*(q*q*(h(m,q)-h(mass,q0))+(mass*mass-sx)*q0*q0*dh(mass,q0));
754 return f;
755}
756double EvtD0ToKpipi0pi0::d(double mass, double q0)const
757{
758 double d = 3.0/pi*mpi*mpi/(q0*q0)*log((mass+2*q0)/(2*mpi))+mass/(2*pi*q0) - (mpi*mpi*mass)/(pi*pow(q0,3));
759 return d;
760}
761EvtComplex EvtD0ToKpipi0pi0::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r, int l)const
762{
763 EvtComplex ci(0,1);
764 EvtComplex prop=1.0/(mass*mass-sa-ci*mass*width*wid(mass,sa,sb,sc,r,l));
765 return prop;
766}
767EvtComplex EvtD0ToKpipi0pi0::propagatorGS(double mass, double width, double sa, double sb, double sc, double r, int l)const
768{
769 EvtComplex ci(0,1);
770 double q = Qabcs(sa,sb,sc);
771 double sa0 = mass*mass;
772 double q0 = Qabcs(sa0,sb,sc);
773 q = sqrt(q);
774 q0 = sqrt(q0);
775 EvtComplex prop = (1+d(mass,q0)*width/mass)/(mass*mass-sa+width*f(mass,sa,q0,q)-ci*mass*width*wid(mass,sa,sb,sc,r,l));
776 return prop;
777}
778double EvtD0ToKpipi0pi0::Flatte_rhoab(double sa, double sb, double sc)const
779{
780 double q = Qabcs(sa,sb,sc);
781 double rho = sqrt(q/sa);
782 return rho;
783}
784EvtComplex EvtD0ToKpipi0pi0::propagatorFlatte(double mass, double width, double sx, double *sb, double *sc)const
785{
786 EvtComplex ci(0,1);
787 double rho1 = Flatte_rhoab(sx,sb[0],sc[0]);
788 double rho2 = Flatte_rhoab(sx,sb[1],sc[1]);
789 EvtComplex prop = 1.0/(mass*mass-sx-ci*(g1*g1*rho1+g2*g2*rho2));
790 return prop;
791}
792double EvtD0ToKpipi0pi0::dot(double *a1, double *a2)const
793{
794 double dot = 0;
795 for(int i=0; i!=4; i++){
796 dot += a1[i]*a2[i]*G[i][i];
797 }
798 return dot;
799}
800double EvtD0ToKpipi0pi0::Qabcs(double sa, double sb, double sc)const
801{
802 double Qabcs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
803 if(Qabcs < 0) Qabcs = 1e-16;
804 return Qabcs;
805}
806double EvtD0ToKpipi0pi0::barrier(double l, double sa, double sb, double sc, double r)const
807{
808 double q = Qabcs(sa,sb,sc);
809 q = sqrt(q);
810 double z = q*r;
811 z = z*z;
812 double F = 1;
813 if(l > 2) F = 0;
814 if(l == 0)
815 F = 1;
816 if(l == 1){
817 F = sqrt((2*z)/(1+z));
818 }
819 if(l == 2){
820 F = sqrt((13*z*z)/(9+3*z+z*z));
821 }
822 return F;
823}
824void EvtD0ToKpipi0pi0::calt1(double daug1[], double daug2[], double t1[]) const
825{
826 double p, pq;
827 double pa[4], qa[4];
828 for(int i=0; i!=4; i++){
829 pa[i] = daug1[i] + daug2[i];
830 qa[i] = daug1[i] - daug2[i];
831 }
832 p = dot(pa,pa);
833 pq = dot(pa,qa);
834 for(int i=0; i!=4; i++){
835 t1[i] = qa[i] - pq/p*pa[i];
836 }
837}
838void EvtD0ToKpipi0pi0::calt2(double daug1[], double daug2[], double t2[][4]) const
839{
840 double p,r;
841 double pa[4], t1[4];
842 calt1(daug1,daug2,t1);
843 r = dot(t1,t1);
844 for(int i=0; i!=4; i++){
845 pa[i] = daug1[i] + daug2[i];
846 }
847 p = dot(pa,pa);
848 for(int i=0; i!=4; i++){
849 for(int j=0; j!=4; j++){
850 t2[i][j] = t1[i]*t1[j] - 1.0/3*r*(G[i][j]-pa[i]*pa[j]/p);
851 }
852 }
853}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
mass
double mass
Definition: CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition: KKsem.h:33
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/Evt3Rank3C.cc:175
t
TTree * t
Definition: binning.cxx:23
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtDecayBase.cc:482
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtParticle.cc:1071