BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKKpipi0.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: EvtDToKKpipi0.cc
11// the necessary file: EvtDToKKpipi0.hh
12//
13// Description: D+ -> K- K+ pi+ pi0
14//
15// Modification history:
16// Liaoyuan Dong, Zhenxuan Li May 26 16:55:36 2024 Module created
17//------------------------------------------------------------------------
18#include "EvtGenBase/EvtPatches.hh"
19#include "EvtGenBase/EvtParticle.hh"
20#include "EvtGenBase/EvtGenKine.hh"
21#include "EvtGenBase/EvtPDL.hh"
22#include "EvtGenBase/EvtReport.hh"
23#include "EvtGenModels/EvtDToKKpipi0.hh"
24#include "EvtGenBase/EvtComplex.hh"
25#include "EvtGenBase/EvtDecayTable.hh"
26#include <stdlib.h>
27#include <stdio.h>
28#include <iostream>
29#include <cmath>
30#include "TSpline.h"
31#include "TComplex.h"
32#include "TMath.h"
33#include "TMatrixD.h"
34using namespace std;
35
36EvtDToKKpipi0::~EvtDToKKpipi0() {}
37
38void EvtDToKKpipi0::getName(std::string& model_name){
39 model_name="DToKKpipi0";
40}
41
45
46void EvtDToKKpipi0::init(){
47
48 checkNArg(0);
49 checkNDaug(4);
50 checkSpinParent(EvtSpinType::SCALAR);
51
52 mD = 1.86966; //D+ mass
53 // mD = 1.86483; //D0 mass
54 mDs = 1.9683;
55 rRes = 3.0;
56 rD = 5.0;
57 metap = 0.95778;
58 mkstr = 0.89594;
59 mk0 = 0.497614;
60 mass_Kaon = 0.49368;
61 mass_Pion = 0.13957;
62 mass_Eta = 0.547862;
63 mass_Etap = 0.95778;
64 mass_Pi0 = 0.1349766;
65 math_pi = 3.1415926;
66 mrho = 0.77511;
67 Grho = 1.2037e-01;
68 mrho1450 = 1.465;
69 Grho1450 = 0.400;
70 ma0 = 0.990;
71 Ga0 = 0.0756;
72 ma01450 = 1.474;
73 Ga01450 = 0.265;
74 meta = 0.547862;
75 ma2 = 1.3177;
76 ma20 = 1.698;
77 Ga2 = 0.1111;
78 Ga20 = 0.265;
79 mf980 = 0.965;
80 Gf980 = 0.0756;
81 meta1405 = 1.4088;
82 Geta1405 = 0.0501;
83 //---------------------------------
84 mphi = 1.019461;
85 Gphi = 0.004249;
86 // Gphi = 6.6912e-03;
87 mnKstr = 0.89555; GnKstr = 0.0473;
88 mcKstr = 0.89167; GcKstr = 0.0514;
89 mK1270 = 1.28981;
90 mK1400 = 1.403;
91 mf1285 = 1.2819; Gf1285 = 0.0227;
92 mf1420 = 1.4263; Gf1420 = 0.0545;
93 gpieta2 = 0.341;
94 gKK2 = 0.892*0.341;
95
96 phi[0] = 0.0; rho[0] = 1.0; //D2VV 12340 phi rho S
97 phi[1] = -2.4183895412; rho[1] = 0.3660444327; //D2VV 12340 phi rho P
98 phi[2] = 3.9126761898; rho[2] = 0.9905707322; //D2VV 12340 phi rho D
99 phi[3] = -3.4581051622; rho[3] = 0.4970069076;
100 phi[4] = 0.6059414627; rho[4] = 0.4744002304;
101 phi[5] = 0.7501566868; rho[5] = 0.6407798868;
102 phi[6] = -1.3470518401; rho[6] = -2.2701523251;
103 phi[7] = -0.9948321319; rho[7] = -1.0841075539;
104 phi[8] = -0.7144319026; rho[8] = 1.3650576749;
105 phi[9] = -0.0298815830; rho[9] = -0.3101468271;
106 phi[10]= -1.0962771498; rho[10]= 0.3799755719;
107 phi[11]= 0.2459805626; rho[11]= -0.3419295418;
108 phi[12]= -1.4937738501; rho[12]= -0.4275950103;
109
110 sp1270 = new TSpline3("sp1270", s, width1270, 1206);
111 sp1400 = new TSpline3("sp1400", s, width1400, 1206);
112
113 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
114 int EE[4][4][4][4] =
115 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
116 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
117 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
118 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
119 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
120 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
121 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
122 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
123 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
124 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
125 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
126 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
127 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
128 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
129 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
130 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
131 for (int i=0; i<4; i++) {
132 for (int j=0; j<4; j++) {
133 G[i][j] = GG[i][j];
134 for (int k=0; k<4; k++) {
135 for (int l=0; l<4; l++) {
136 E[i][j][k][l] = EE[i][j][k][l];
137 }
138 }
139 }
140 }
141}
142
143void EvtDToKKpipi0::initProbMax() {
144 setProbMax(301000);
145}
146
147void EvtDToKKpipi0::decay( EvtParticle *p){
148 //-----------for max value------------------
149 /*double maxprob = 0.0;
150 for(int ir=0;ir<=60000000;ir++){
151 p->initializePhaseSpace(getNDaug(),getDaugs());
152 EvtVector4R p1 = p->getDaug(0)->getP4();
153 EvtVector4R p2 = p->getDaug(1)->getP4();
154 EvtVector4R p3 = p->getDaug(2)->getP4();
155 EvtVector4R p4 = p->getDaug(3)->getP4();
156 double value;
157 double P1[4],P2[4],P3[4],P4[4];
158 mother_c=EvtPDL::getStdHep(p->getId());
159 if(mother_c==411){
160 P1[0] = p1.get(0); P2[0] = p2.get(0); P3[0] = p3.get(0); P4[0] = p4.get(0);
161 P1[1] = p1.get(1); P2[1] = p2.get(1); P3[1] = p3.get(1); P4[1] = p4.get(1);
162 P1[2] = p1.get(2); P2[2] = p2.get(2); P3[2] = p3.get(2); P4[2] = p4.get(2);
163 P1[3] = p1.get(3); P2[3] = p2.get(3); P3[3] = p3.get(3); P4[3] = p4.get(3);
164 }else if(mother_c==-411){
165 P1[0] = p1.get(0); P2[0] = p2.get(0); P3[0] = p3.get(0); P4[0] = p4.get(0);
166 P1[1] = -1.0*p1.get(1); P2[1] = -1.0*p2.get(1); P3[1] = -1.0*p3.get(1); P4[1] = -1.0*p4.get(1);
167 P1[2] = -1.0*p1.get(2); P2[2] = -1.0*p2.get(2); P3[2] = -1.0*p3.get(2); P4[2] = -1.0*p4.get(2);
168 P1[3] = -1.0*p1.get(3); P2[3] = -1.0*p2.get(3); P3[3] = -1.0*p3.get(3); P4[3] = -1.0*p4.get(3);
169 }
170 calPDF(P1, P2, P3, P4, value);
171 if(value>maxprob) {
172 maxprob=value;
173 std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
174 }
175 }
176 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
177 return;*/
178 //-----------------------------------------------
179 p->initializePhaseSpace(getNDaug(),getDaugs());
180 EvtVector4R p1 = p->getDaug(0)->getP4();
181 EvtVector4R p2 = p->getDaug(1)->getP4();
182 EvtVector4R p3 = p->getDaug(2)->getP4();
183 EvtVector4R p4 = p->getDaug(3)->getP4();
184 mother_c=EvtPDL::getStdHep(p->getId());
185
186 double P1[4],P2[4],P3[4],P4[4];
187 if(mother_c==411){
188 P1[0] = p1.get(0); P2[0] = p2.get(0); P3[0] = p3.get(0); P4[0] = p4.get(0);
189 P1[1] = p1.get(1); P2[1] = p2.get(1); P3[1] = p3.get(1); P4[1] = p4.get(1);
190 P1[2] = p1.get(2); P2[2] = p2.get(2); P3[2] = p3.get(2); P4[2] = p4.get(2);
191 P1[3] = p1.get(3); P2[3] = p2.get(3); P3[3] = p3.get(3); P4[3] = p4.get(3);
192 }else if(mother_c==-411){
193 P1[0] = p1.get(0); P2[0] = p2.get(0); P3[0] = p3.get(0); P4[0] = p4.get(0);
194 P1[1] = -1.0*p1.get(1); P2[1] = -1.0*p2.get(1); P3[1] = -1.0*p3.get(1); P4[1] = -1.0*p4.get(1);
195 P1[2] = -1.0*p1.get(2); P2[2] = -1.0*p2.get(2); P3[2] = -1.0*p3.get(2); P4[2] = -1.0*p4.get(2);
196 P1[3] = -1.0*p1.get(3); P2[3] = -1.0*p2.get(3); P3[3] = -1.0*p3.get(3); P4[3] = -1.0*p4.get(3);
197 }
198
199 double value;
200 calPDF(P1, P2, P3, P4, value);
201 setProb(value);
202 return;
203}
204
205double EvtDToKKpipi0::calPDF(double P1[], double P2[], double P3[], double P4[], double & Result) {
206
207 double s12, s13, s14, s23, s24, s34;
208 double s123, s124, s134, s234;
209
210 double P12[4], P13[4], P14[4], P23[4], P24[4], P34[4];
211 double P123[4], P124[4], P134[4], P234[4];
212 for(int i = 0; i < 4; i++){
213 P12[i] = P1[i] + P2[i];
214 P13[i] = P1[i] + P3[i];
215 P14[i] = P1[i] + P4[i];
216 P23[i] = P2[i] + P3[i];
217 P24[i] = P2[i] + P4[i];
218 P34[i] = P3[i] + P4[i];
219 P123[i] = P12[i] + P3[i];
220 P124[i] = P12[i] + P4[i];
221 P134[i] = P13[i] + P4[i];
222 P234[i] = P23[i] + P4[i];
223 }
224
225 s12 = dot(P12,P12);
226 s13 = dot(P13,P13);
227 s14 = dot(P14,P14);
228 s23 = dot(P23,P23);
229 s24 = dot(P24,P24);
230 s34 = dot(P34,P34);
231
232 s123 = dot(P123,P123);
233 s124 = dot(P124,P124);
234 s134 = dot(P134,P134);
235 s234 = dot(P234,P234);
236 TComplex PDF[100], cof, pdf, module, pro1, pro2, cpro1, cpro2;
237 //---------
238 double wid1270, wid1400;
239 wid1270 = sp1270->Eval(s134);
240 wid1400 = sp1400->Eval(s134);
241 //---------
242
243 pro1 = propagatorRBW(mphi,Gphi,s12,sck,sck,3.0,1,1);
244 pro2 = propagatorGS(mrho,Grho,s34,scpi,snpi,3.0,1);
245 PDF[0] = Spin_factor(P1, P2, P3, P4, 0, mphi, mrho, 3.0, 5.0, 0)*pro1*pro2; // phi rho S-wave
246 PDF[1] = Spin_factor(P1, P2, P3, P4, 1, mphi, mrho, 3.0, 5.0, 0)*pro1*pro2; // phi rho P-wave
247 PDF[2] = Spin_factor(P1, P2, P3, P4, 2, mphi, mrho, 3.0, 5.0, 0)*pro1*pro2; // phi rho D-wave
248
249 pro1 = propagatorRBW(mnKstr,GnKstr,s13,sck,scpi,3.0,1,1);
250 pro2 = propagatorRBW(mcKstr,GcKstr,s24,sck,snpi,3.0,1,1);
251 PDF[3] = Spin_factor(P1, P3, P2, P4, 0, mnKstr, mcKstr, 3.0, 5.0, 0)*pro1*pro2; // K*0bar K*+ S-wave
252 PDF[4] = Spin_factor(P1, P3, P2, P4, 1, mnKstr, mcKstr, 3.0, 5.0, 0)*pro1*pro2; // K*0bar K*+ P-wave
253 PDF[5] = Spin_factor(P1, P3, P2, P4, 2, mnKstr, mcKstr, 3.0, 5.0, 0)*pro1*pro2; // K*0bar K*+ D-wave
254
255 pro1 = propagatorGS(mrho,Grho,s34,scpi,snpi,3.0,1);
256 pro2 = propagatorRBW(mK1270,wid1270,s134,s34,sck,3.0,1,0);
257 PDF[6] = Spin_factor(P3, P4, P1, P2, 0, mrho, mK1270, 3.0, 5.0, 1)*pro1*pro2; // K1270 pi, K1270->rho K S-wave
258
259 pro1 = propagatorRBW(mnKstr,GnKstr,s13,sck,scpi,3.0,1,1);
260 pro2 = propagatorRBW(mK1400,wid1400,s134,s13,snpi,3.0,1,0);
261 cpro1 = propagatorRBW(mcKstr,GcKstr,s14,sck,snpi,3.0,1,1);
262 cpro2 = propagatorRBW(mK1400,wid1400,s134,s14,scpi,3.0,1,0);
263 PDF[7] = Spin_factor(P1, P3, P4, P2, 0, mnKstr, mK1400, 3.0, 5.0, 1)*pro1*pro2 + Spin_factor(P1, P4, P3, P2, 0, mcKstr, mK1400, 3.0, 5.0, 1)*cpro1*cpro2; //K1400 pi, K1400->Kstr pi
264
265 pro1 = propagatorGS(mrho,Grho,s34,scpi,snpi,3.0,1);
266 double sb[2], sc[2], pro_tmp[2];
267 //sb[0] = snpi; sc[0] = seta;
268 //sb[1] = sck; sc[1] = sck;
269 //pro2 = propagatorFlatte(ma0, Ga0, s12, sb, sc, gpieta2, gKK2);
270 propagatorFlatte(ma0, Ga0, s12, snpi, seta, 0,pro_tmp);
271 pro2 = TComplex(pro_tmp[0],pro_tmp[1]);
272 PDF[8] = Spin_factor(P3, P4, P1, P2, 0, mrho, ma0, 3.0, 5.0, 10)*pro1*pro2;// rho a0
273
274 pro1 = propagatorRBW(mcKstr,GcKstr,s14,sck,snpi,3.0,1,1);
275 pro2 = propagatorRBW(mf1420,Gf1420,s124,s14,sck,3.0,0,1);
276 cpro1 = propagatorRBW(mcKstr,GcKstr,s24,sck,snpi,3.0,1,1);
277 cpro2 = propagatorRBW(mf1420,Gf1420,s124,s24,sck,3.0,0,1);
278 PDF[9] = Spin_factor(P1, P4, P2, P3, 0, mcKstr, mf1420, 3.0, 5.0, 1)*pro1*pro2 + Spin_factor(P2, P4, P1, P3, 0, mcKstr, mf1420, 3.0, 5.0, 1)*cpro1*cpro2; //f1420 -> Kstr K
279
280 pro1 = propagatorRBW(mcKstr,GcKstr,s14,sck,snpi,3.0,1,1);
281 pro2 = propagatorRBW(mf1285,Gf1285,s124,s14,sck,3.0,0,1);
282 cpro1 = propagatorRBW(mcKstr,GcKstr,s24,sck,snpi,3.0,1,1);
283 cpro2 = propagatorRBW(mf1285,Gf1285,s124,s24,sck,3.0,0,1);
284 PDF[10] = Spin_factor(P1, P4, P2, P3, 0, mcKstr, mf1285, 3.0, 5.0, 1)*pro1*pro2 + Spin_factor(P2, P4, P1, P3, 0, mcKstr, mf1285, 3.0, 5.0, 1)*cpro1*cpro2; //f1285 -> Kstr K
285
286 //pro1 = propagatorFlatte(ma0, Ga0, s12, sb, sc, gpieta2, gKK2);
287 propagatorFlatte(ma0, Ga0, s12, snpi, seta, 0,pro_tmp);
288 pro1 = TComplex(pro_tmp[0],pro_tmp[1]);
289 pro2 = propagatorRBW(meta1405,Geta1405,s124,s12,snpi,3.0,1,1);
290 PDF[11] = pro1*pro2; //eta1405->a0 pi
291
292 pro1 = propagatorRBW(mcKstr,GcKstr,s14,sck,snpi,3.0,1,1);
293 pro2 = TComplex(pro_tmp[0],pro_tmp[1]);
294 pro2 = propagatorRBW(meta1405,Geta1405,s124,s14,sck,3.0,1,1);
295 cpro1 = propagatorRBW(mcKstr,GcKstr,s24,sck,snpi,3.0,1,1);
296 cpro2 = propagatorRBW(meta1405,Geta1405,s124,s24,sck,3.0,1,1);
297 PDF[12] = Spin_factor(P1, P4, P2, P3, 1, mcKstr, meta1405, 3.0, 5.0, 21)*pro1*pro2 + Spin_factor(P2, P4, P1, P3, 1, mcKstr, meta1405, 3.0, 5.0, 21)*cpro1*cpro2; //eta1405 -> Kstr K
298
299 pdf = TComplex(0,0);
300 for(Int_t i=0; i<13; i++){
301 cof = TComplex(rho[i]*TMath::Cos(phi[i]),rho[i]*TMath::Sin(phi[i]));
302 //cout<<i<<" "<<(cof*PDF[i]).Re()<<" "<<(cof*PDF[i]).Im()<<endl;
303 pdf += cof*PDF[i];
304 }
305 module = pdf*TComplex::Conjugate(pdf);
306 double value;
307 value = module.Re();
308 // cout<<value<<endl;
309 // return (value <= 0) ? 1e-20 : value;
310 Result = (value <= 0) ? 1e-20 : value;
311}
312
313TComplex EvtDToKKpipi0::propogator(double mass, double width, double sx)const
314{
315 TComplex ci(0,1);
316 TComplex prop=1.0/(mass*mass-sx-ci*mass*width);
317 return prop;
318}
319double EvtDToKKpipi0::wid(double mass, double sa, double sb, double sc, double r, int l)const
320{
321 double widm(0.), q(0.), q0(0.);
322 double sa0 = mass*mass;
323 double m = sqrt(sa);
324 q = Qabcs(sa,sb,sc);
325 q0 = Qabcs(sa0,sb,sc);
326 double z,z0;
327 z = q*r*r;
328 z0 = q0*r*r;
329 double t = q/q0;
330 double F(0.);
331 if(l == 0) F = 1;
332 if(l == 1) F = sqrt((1+z0)/(1+z));
333 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
334 widm = pow(t,l+0.5)*mass/m*F*F;
335 return widm;
336}
337double EvtDToKKpipi0::h(double m, double q)const
338{
339 double h(0.);
340 h = 2/math_pi*q/m*log((m+2*q)/(0.13957 + 0.134976));
341 return h;
342}
343double EvtDToKKpipi0::dh(double mass, double q0)const
344{
345 double dh = h(mass,q0)*(1.0/(8*q0*q0)-1.0/(2*mass*mass))+1.0/(2*math_pi*mass*mass);
346 return dh;
347}
348double EvtDToKKpipi0::f(double mass, double sx, double q0, double q)const
349{
350 double m = sqrt(sx);
351 double f = mass*mass/(pow(q0,3))*(q*q*(h(m,q)-h(mass,q0))+(mass*mass-sx)*q0*q0*dh(mass,q0));
352 return f;
353}
354double EvtDToKKpipi0::d(double mass, double q0)const
355{
356 double cmpi = 0.5*(0.13957 + 0.134976);
357 // double cmpi = 0.13957;
358 double d = 3.0/math_pi*cmpi*cmpi/(q0*q0)*log((mass+2*q0)/(2*cmpi))+mass/(2*math_pi*q0)
359 -(cmpi*cmpi*mass)/(math_pi*pow(q0,3));
360 return d;
361}
362TComplex EvtDToKKpipi0::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r, int l, int isvaried)const
363{
364 TComplex ci(0,1);
365 TComplex prop;
366 if(isvaried == 1) prop=1.0/(mass*mass-sa-ci*mass*width*wid(mass,sa,sb,sc,r,l));
367 if(isvaried == 0) prop=1.0/(mass*mass-sa-ci*mass*width);
368 return prop;
369}
370TComplex EvtDToKKpipi0::propagatorGS(double mass, double width, double sa, double sb, double sc, double r, int l)const
371{
372 TComplex ci(0,1);
373 double q = Qabcs(sa,sb,sc);
374 double sa0 = mass*mass;
375 double q0 = Qabcs(sa0,sb,sc);
376 q = sqrt(q);
377 q0 = sqrt(q0);
378 TComplex 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));
379 // TComplex prop = (1+d(mass,q0)*width/mass)/(mass*mass-sa+width*f(mass,sa,q0,q)-ci*mass*width);
380 return prop;
381}
382TComplex EvtDToKKpipi0::Flatte_rhoab(double sa, double sb, double sc)const
383{
384 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
385 TComplex rho;
386 TComplex ci(0,1);
387 if(q>0){
388 rho = 2.0*sqrt(q/sa)*(TComplex::One());
389 }
390 if(q<0){
391 rho = 2.0*sqrt(-q/sa)*ci;
392 }
393 return rho;
394}
395void EvtDToKKpipi0::Com_Divide(double a1[2], double a2[2], double res[2]){
396 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
397 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
398}
399void EvtDToKKpipi0::propagatorFlatte(double mass, double width, double sa, double sb, double sc, int r, double prop[2])
400{
401 //r == 0, neutral a0
402 //r == 1, charged a0
403 double q, qKsK;
404 double rhoab[2], rhoKsK[2];
405 q = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
406 if(r == 0) qKsK = 0.25*sa - 0.49368*0.49368;
407 if(r == 1) qKsK = 0.25*(sa + 3.899750596e-03)*(sa + 3.899750596e-03)/sa - 0.497614*0.497614;
408 if(q>0){
409 rhoab[0] = 2*sqrt(q/sa);
410 rhoab[1] = 0;
411 }
412 if(q<0){
413 rhoab[0] = 0;
414 rhoab[1] = 2*sqrt(-q/sa);
415 }
416 if(qKsK>0){
417 rhoKsK[0] = 2*sqrt(qKsK/sa);
418 rhoKsK[1] = 0;
419 }
420 if(qKsK<0){
421 rhoKsK[0] = 0;
422 rhoKsK[1] = 2*sqrt(-qKsK/sa);
423 }
424 double a[2], b[2];
425 a[0] = 1;
426 a[1] = 0;
427 b[0] = mass*mass - sa + 0.341*rhoab[1] + 0.892*0.341*rhoKsK[1];
428 b[1] = - (0.341*rhoab[0] + 0.892*0.341*rhoKsK[0]);
429 Com_Divide(a,b,prop);
430}
431double EvtDToKKpipi0::Spin_factor(double P1[], double P2[], double P3[], double P4[], Int_t Ang, double mass1, double mass2, double rRes, double rD, Int_t mode)
432{
433 //mode == 0, D->VV
434 //mode == 1, D->AP, A->VP, V->PP
435 //mode == 2, D->VP, V->VP, V->PP
436
437 //mode == 10, D->VS, V->PP, S->PP
438 //mode == 11, D->TS, T->PP, S->PP
439 //mode == 12, PHSP
440 //mode == 20, D->AP, A->SP, S->PP
441 //mode == 21, D->P'P, P'->VP, V->PP
442
443 double pR1[4], pR2[4], pD[4];
444 double temp_PDF, v_re;
445 temp_PDF = 0.0;
446 v_re = 0.0;
447 double mD = 1.86966;
448 double B[3], s1, s2, s3, s4, sR1, sR2, sD;
449 sD = mD*mD;
450 double t1D[4], t2D[4][4];
451 s1 = dot(P1,P1);
452 s2 = dot(P2,P2);
453 s3 = dot(P3,P3);
454 s4 = dot(P4,P4);
455 double tV1[4], tV2[4];
456
457 //=============================
458 if(mode == 0){
459 //VV
460 for(int i=0; i<4; i++){
461 pR1[i] = P1[i] + P2[i];
462 pR2[i] = P3[i] + P4[i];
463 pD[i] = pR1[i] + pR2[i];
464 }
465 sR1 = dot(pR1,pR1);
466 sR2 = dot(pR2,pR2);
467 calt1(P1,P2,tV1);
468 calt1(P3,P4,tV2);
469 calt1(pR1,pR2,t1D);
470 calt2(pR1,pR2,t2D);
471 B[0] = barrier(1,sR1,s1,s2,rRes,mass1);
472 B[1] = barrier(1,sR2,s3,s4,rRes,mass2);
473
474 if(Ang == 0){
475 B[2] = 1.0;
476 // for(int i=0; i<4; i++){
477 // temp_PDF += tV1[i]*tV2[i]*G[i][i];
478 // }
479 for(int i=0; i<4; i++){
480 for(int j=0; j<4; j++){
481 temp_PDF += (-G[i][j] + pD[i]*pD[j]/sD)*tV1[i]*tV2[j]*G[i][i]*G[j][j];
482 }
483 }
484 }
485 if(Ang == 1){
486 B[2] = barrier(1,sD,sR1,sR2,rD,mD);
487 for(int i=0; i<4; i++){
488 for(int j=0; j<4; j++){
489 for(int k=0; k<4; k++){
490 for(int l=0; l<4; l++){
491 temp_PDF += E[i][j][k][l]*pD[i]*t1D[j]*tV1[k]*tV2[l]*G[i][i]*G[j][j]*G[k][k]*G[l][l];
492 }
493 }
494 }
495 }
496 }
497 if(Ang == 2){
498 B[2] = barrier(2,sD,sR1,sR2,rD,mD);
499 for(int i=0; i<4; i++){
500 for(int j=0; j<4; j++){
501 // temp_PDF += pD[i]*pD[j]*tV1[i]*tV2[j]*G[i][i]*G[j][j];
502 temp_PDF += t2D[i][j]*tV1[i]*tV2[j]*G[i][i]*G[j][j];
503 }
504 }
505 }
506 }
507 //=============================
508 if(mode == 1){
509 // AP
510 for(int i=0; i<4; i++){
511 pR1[i] = P1[i] + P2[i];
512 pR2[i] = pR1[i] + P3[i];
513 pD[i] = pR2[i] + P4[i];
514 }
515 double tV[4];
516 sR1 = dot(pR1,pR1);
517 sR2 = dot(pR2,pR2);
518 calt1(P1,P2,tV);
519 calt1(pR2,P4,t1D);
520 B[0] = barrier(1,sR1,s1,s2,rRes,mass1);
521 B[2] = barrier(1,sD,sR2,s4,rD,mD);
522 if(Ang == 0){
523 B[1] = 1.0;
524 for(int i=0; i<4; i++){
525 for(int j=0; j<4; j++){
526 temp_PDF += t1D[i]*(-G[i][j] + pR2[i]*pR2[j]/sR2)*tV[j]*G[i][i]*G[j][j];
527 }
528 }
529 }
530 if(Ang == 2){
531 B[1] = barrier(2,sR2,sR1,s3,rRes,mass2);
532 double t2A[4][4];
533 calt2(pR1,P3,t2A);
534 for(int i=0; i<4; i++){
535 for(int j=0; j<4; j++){
536 temp_PDF += t1D[i]*t2A[i][j]*tV[j]*G[i][i]*G[j][j];
537 }
538 }
539 }
540 }
541 //=============================
542 if(mode == 2){
543 //VP
544 for(int i=0; i<4; i++){
545 pR1[i] = P1[i] + P2[i];
546 pR2[i] = pR1[i] + P3[i];
547 pD[i] = pR2[i] + P4[i];
548 }
549 sR1 = dot(pR1,pR1);
550 sR2 = dot(pR2,pR2);
551 B[0] = barrier(1,sR1,s1,s2,rRes,mass1);
552 B[1] = barrier(1,sR2,sR1,s3,rRes,mass2);
553 B[2] = barrier(1,sD,sR2,s4,rD,mD);
554 for(int i=0; i<4; i++){
555 for(int j=0; j<4; j++){
556 for(int k=0; k<4; k++){
557 for(int l=0; l<4; l++){
558 temp_PDF += E[i][j][k][l]*pR2[i]*(pR1[j] - P3[j])*P4[k]*(P1[l] - P2[l])*G[i][i]*G[j][j]*G[k][k]*G[l][l];
559 }
560 }
561 }
562 }
563 }
564 //=============================
565 if(mode == 10){
566 //VS
567 for(int i=0; i<4; i++){
568 pR1[i] = P1[i] + P2[i];
569 pR2[i] = P3[i] + P4[i];
570 pD[i] = pR1[i] + pR2[i];
571 }
572 sR1 = dot(pR1,pR1);
573 sR2 = dot(pR2,pR2);
574 calt1(P1,P2,tV1);
575 calt1(pR1,pR2,t1D);
576 B[0] = barrier(1,sR1,s1,s2,rRes,mass1);
577 B[1] = 1.0;
578 B[2] = barrier(1,sD,sR1,sR2,rD,mD);
579 for(int i=0; i<4; i++){
580 temp_PDF += t1D[i]*tV1[i]*G[i][i];
581 }
582
583 }
584 //=============================
585 if(mode == 11){
586 // TS
587 for(int i=0; i<4; i++){
588 pR1[i] = P1[i] + P2[i];
589 pR2[i] = P3[i] + P4[i];
590 pD[i] = pR1[i] + pR2[i];
591 }
592 sR1 = dot(pR1,pR1);
593 sR2 = dot(pR2,pR2);
594 double tT[4][4], t2D[4][4];
595 calt2(P1,P2,tT);
596 calt2(pR1,pR2,t2D);
597 B[0] = barrier(2,sR1,s1,s2,rRes,mass1);
598 B[1] = 1.0;
599 B[2] = barrier(2,sD,sR1,sR2,rD,mD);
600 for(int i=0; i<4; i++){
601 for(int j=0; j<4; j++){
602 temp_PDF += t2D[i][j]*tT[i][j]*G[i][i]*G[j][j];
603 }
604 }
605 }
606 //=============================
607 if(mode == 20){
608 //AP, A->SP
609 for(int i=0; i<4; i++){
610 pR1[i] = P1[i] + P2[i];
611 pR2[i] = pR1[i] + P3[i];
612 pD[i] = pR2[i] + P4[i];
613 }
614 double tA[4];
615 sR1 = dot(pR1,pR1);
616 sR2 = dot(pR2,pR2);
617 calt1(pR1,P3,tA);
618 calt1(pR2,P4,t1D);
619 B[0] = 1.0;
620 B[1] = barrier(1,sR2,sR1,s3,rRes,mass2);
621 B[2] = barrier(1,sD,sR2,s4,rD,mD);
622 for(int i=0; i<4; i++){
623 temp_PDF += tA[i]*t1D[i]*G[i][i];
624 }
625 }
626 //=============================
627 if(mode == 21){
628 //P' P, P'->VP, V->PP
629 for(int i=0; i<4; i++){
630 pR1[i] = P1[i] + P2[i];
631 pR2[i] = pR1[i] + P3[i];
632 pD[i] = pR2[i] + P4[i];
633 }
634 double tVP[4];
635 sR1 = dot(pR1,pR1);
636 sR2 = dot(pR2,pR2);
637 calt1(P1,P2,tVP);
638 B[0] = barrier(1,sR1,s1,s2,rRes,mass1);
639 B[1] = barrier(1,sR2,sR1,s3,rRes,mass2);
640 B[2] = 1.0;
641 for(int i=0; i<4; i++){
642 temp_PDF += P3[i]*tVP[i]*G[i][i];
643 }
644 }
645 //=============================
646 v_re = temp_PDF*B[0]*B[1]*B[2];
647 return v_re;
648}
649double EvtDToKKpipi0::CalRho4pi(double s){
650
651 if(s>=1.){
652 return sqrt((s-16.*mass_Pion*mass_Pion)/s);
653 }
654 else{
655 double gam = 0.0005;
656 double s2 = s*s;
657 double s3 = s2*s;
658 double s4 = s2*s2;
659 double s5 = s4*s;
660 double s6 = s5*s;
661
662 gam -= 0.0193*s;
663 gam += 0.1385*s2;
664 gam -= 0.2084*s3;
665 gam -= 0.2974*s4;
666 gam += 0.1366*s5;
667 gam += 1.0789*s6;
668
669 return gam;
670 }
671}
672
673TComplex EvtDToKKpipi0::ResonanceSkm(double m2){
674
675 //all par value
676 double g11 = 0.22889, g12 = -0.55377, g13 = 0.00, g14 = -0.39899, g15 = -0.34639;
677 double g21 = 0.94128, g22 = 0.55095, g23 = 0.00, g24 = 0.39065, g25 = 0.31503;
678 double g31 = 0.36856, g32 = 0.23888, g33 = 0.55639, g34 = 0.18340, g35 = 0.18681;
679 double g41 = 0.33650, g42 = 0.40907, g43 = 0.85679, g44 = 0.19906, g45 = -0.00984;
680 double g51 = 0.18171, g52 = -0.17558, g53 = -0.79658, g54 = -0.00355, g55 = 0.22358;
681
682 double fs11 = 0.23399, fs12 = 0.15044, fs13 =-0.20545, fs14 = 0.32825, fs15 = 0.35412;
683 double m_1 = 0.65100, m_2 = 1.20360, m_3 = 1.55817, m_4 = 1.21000, m_5 = 1.82206;
684 double ss0 = -3.92637;
685 double sp0 = -3.0;//v0
686 //double sp0 = -0.07;//v1
687 double sA0 = -0.15;
688 double sA = 1.0;
689
690 double 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);
691 double 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);
692 double 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);
693 double 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);
694 double 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);
695
696 double km21 = km12, km31 = km13, km41 = km14, km51 = km15;
697
698 double 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);
699 double 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);
700 double 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);
701
702 double km32 = km23, km42 = km24, km52 = km25;
703 double 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);
704 double 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);
705 double km43 = km34, km53 = km35;
706 double 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);
707 double km54 = km45;
708
709 double 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);
710 double 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);
711 double 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);
712 double 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);
713 TComplex fp11(16.55, -3.87); TComplex beta1(8.89, 3.36);
714 TComplex fp12(-13.99, 5.12); TComplex beta2(16.20, 5.77);
715 TComplex fp13(-77.34, -80.04); TComplex beta3(-21.60, 27.41);
716 TComplex fp14(-28.52, -3.20); TComplex beta4(-40.15, 35.36);
717 TComplex fp15(30.47, 5.70); TComplex beta5(30.0, -43.09);
718
719 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);
720 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);
721 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);
722 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);
723 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);
724
725 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);
726 MI.UnitMatrix();
727 RhoA.UnitMatrix();
728 RhoB.UnitMatrix();
729
730 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;
731 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;
732 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;
733 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;
734 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;
735
736 TMatrixDRow(RhoA,0)(0) = sqrt(1.0-4.0*mass_Pion*mass_Pion/m2);
737 TMatrixDRow(RhoB,0)(0) = 0.0;
738 if ( (1.0-4.0*mass_Kaon*mass_Kaon/m2) > 0) {
739 TMatrixDRow(RhoA,1)(1) = sqrt(1.0-4.0*mass_Kaon*mass_Kaon/m2);
740 TMatrixDRow(RhoB,1)(1) = 0.0;
741 }
742 else {
743 TMatrixDRow(RhoA,1)(1) = 0.0;
744 TMatrixDRow(RhoB,1)(1) = sqrt(4.0*mass_Kaon*mass_Kaon/m2-1.0);
745 }
746
747 TMatrixDRow(RhoA,2)(2) = CalRho4pi(m2);
748 TMatrixDRow(RhoB,2)(2) = 0.0;
749 if ( (1.0-4.0*mass_Eta*mass_Eta/m2) > 0) {
750 TMatrixDRow(RhoA,3)(3) = sqrt(1.0-4.0*mass_Eta*mass_Eta/m2);
751 TMatrixDRow(RhoB,3)(3) = 0.0;
752 }else {
753 TMatrixDRow(RhoA,3)(3) = 0.0;
754 TMatrixDRow(RhoB,3)(3) = sqrt(4.0*mass_Eta*mass_Eta/m2-1.0);
755 }
756 TMatrixDRow(RhoA,4)(4) = 0.0;
757 TMatrixDRow(RhoB,4)(4) = sqrt((mass_Eta + mass_Etap)*(mass_Eta + mass_Etap)/m2-1.0);
758
759 MA = MI + KM*RhoB;
760 MB = -1.0*KM*RhoA;
761
762 MA_invt = MA;
763 MA_invt.Invert();
764
765 MRe = MA+MB*MA_invt*MB;
766 MRe.Invert();
767 MIm = MA_invt*MB*MRe;
768
769 TComplex ciR(1.0, 0.0);
770 TComplex ciM(0.0, 1.0);
771 TComplex amp;
772
773 amp = (ciR*TMatrixDRow(MRe,0)(0)-ciM*TMatrixDRow(MIm,0)(0))*P1+
774 (ciR*TMatrixDRow(MRe,0)(1)-ciM*TMatrixDRow(MIm,0)(1))*P2+
775 (ciR*TMatrixDRow(MRe,0)(2)-ciM*TMatrixDRow(MIm,0)(2))*P3+
776 (ciR*TMatrixDRow(MRe,0)(3)-ciM*TMatrixDRow(MIm,0)(3))*P4+
777 (ciR*TMatrixDRow(MRe,0)(4)-ciM*TMatrixDRow(MIm,0)(4))*P5;
778
779 return amp;
780}
781
782double EvtDToKKpipi0::dot(double *a1, double *a2)const
783{
784 double dot = 0;
785 for(Int_t i=0; i!=4; i++){
786 dot += a1[i]*a2[i]*G[i][i];
787 }
788 return dot;
789}
790double EvtDToKKpipi0::Qabcs(double sa, double sb, double sc)const
791{
792 Double_t Qabcs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
793 if(Qabcs < 0) Qabcs = 1e-16;
794 return Qabcs;
795}
796double EvtDToKKpipi0::barrier(double l, double sa, double sb, double sc, double r, double mass)const
797{
798 double sa0 = mass*mass;
799 double q0 = Qabcs(sa0,sb,sc);
800 double z0 = q0*r*r;
801 double q = Qabcs(sa,sb,sc);
802 q = sqrt(q);
803 double z = q*r;
804 z = z*z;
805 double F = 1;
806 if(l > 2) F = 0;
807 if(l == 0)
808 F = 1;
809 if(l == 1){
810 F = sqrt((1+z0)/(1+z));
811 }
812 if(l == 2){
813 F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
814 }
815 return F;
816}
817void EvtDToKKpipi0::calt1(double daug1[], double daug2[], double t1[]) const
818{
819 double p, pq;
820 double pa[4], qa[4];
821 for(int i=0; i!=4; i++){
822 pa[i] = daug1[i] + daug2[i];
823 qa[i] = daug1[i] - daug2[i];
824 }
825 p = dot(pa,pa);
826 pq = dot(pa,qa);
827 for(int i=0; i!=4; i++){
828 t1[i] = qa[i] - pq/p*pa[i];
829 }
830}
831void EvtDToKKpipi0::calt2(Double_t daug1[], Double_t daug2[], Double_t t2[][4]) const
832{
833 double p,r;
834 double pa[4], t1[4];
835 calt1(daug1,daug2,t1);
836 r = dot(t1,t1);
837 for(int i=0; i!=4; i++){
838 pa[i] = daug1[i] + daug2[i];
839 }
840 p = dot(pa,pa);
841 for(int i=0; i!=4; i++){
842 for(int j=0; j!=4; j++){
843 t2[i][j] = t1[i]*t1[j] - 1.0/3*r*(G[i][j]-pa[i]*pa[j]/p);
844 }
845 }
846}
p1
double p1[4]
Definition BabayagaNLO.cxx:43
p2
double p2[4]
Definition BabayagaNLO.cxx:43
mass
double mass
Definition CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
TComplex
complex< double > TComplex
Definition Dalitz.h:14
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
EvtDToKKpipi0::~EvtDToKKpipi0
virtual ~EvtDToKKpipi0()
Definition EvtDToKKpipi0.cc:36
EvtDToKKpipi0::clone
EvtDecayBase * clone()
Definition EvtDToKKpipi0.cc:42
EvtDToKKpipi0::getName
void getName(std::string &name)
Definition EvtDToKKpipi0.cc:38
EvtDToKKpipi0::decay
void decay(EvtParticle *p)
Definition EvtDToKKpipi0.cc:147
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
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition EvtPDL.hh:56
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:112
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
p4
double double double * p4
Definition qcdloop1.h:77
s12
double double double double * s12
Definition qcdloop1.h:77
s23
double double double double double * s23
Definition qcdloop1.h:77
p3
double double * p3
Definition qcdloop1.h:76
m2
double double * m2
Definition qcdloop1.h:75
b
const double b
Definition slope.cxx:9