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