BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKKmunu.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of the EvtGen package developed jointly
4// for the BaBar and CLEO collaborations. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/COPYRIGHT
8// Copyright (C) 1998 Caltech, UCSB
9//
10// Module: EvtDsToKKmunu.cc
11//
12// Description: Routine to handle three-body decays of Ds+/Ds-
13//
14// Modification history:
15// Liaoyuan Dong Feb 12 2024 Module created
16//------------------------------------------------------------------------
17#include "EvtGenBase/EvtPatches.hh"
18#include "EvtGenBase/EvtParticle.hh"
19#include "EvtGenBase/EvtGenKine.hh"
20#include "EvtGenBase/EvtPDL.hh"
21#include "EvtGenBase/EvtReport.hh"
22#include "EvtGenModels/EvtDsToKKmunu.hh"
23#include "EvtGenBase/EvtComplex.hh"
24#include "EvtGenBase/EvtDecayTable.hh"
25#include <stdlib.h>
26
27EvtDsToKKmunu::~EvtDsToKKmunu() {}
28
29void EvtDsToKKmunu::getName(std::string& model_name){
30 model_name="DsToKKmunu";
31}
32
36
37void EvtDsToKKmunu::init(){
38 checkNArg(0);
39 checkNDaug(4);
40 checkSpinParent(EvtSpinType::SCALAR);
41 checkSpinDaughter(0,EvtSpinType::SCALAR);
42 checkSpinDaughter(1,EvtSpinType::SCALAR);
43
44 // std::cout << "EvtDsToKKmunu ==> Initialization !" << std::endl;
45nAmps = 2;
46
47// rS = -11.27; // S-wave
48// rS1 = 0.08;
49// a_delta = 1.58;
50// b_delta = -0.81;
51// m0_1430_S = 1.425;
52// width0_1430_S = 0.270;
53// type[0] = 0;
54
55 mV = 2.1;
56 mA = 2.5;
57 V_0 = 1.58;
58 A1_0 = 1;
59 A2_0 = 0.71;
60
61 m0 = 1.01946; // P-wave phi
62 width0 = 0.004249;
63 rBW = 3.0;
64 rho = 1;
65 phi = 0;
66 type[1] = 1;
67
68 m0_1410 = 1.414; // P-wave K*(1410)
69 width0_1410 = 0.232;
70 rho_1410 = 0.1;
71 phi_1410 = 0.;
72 type[2] = 2;
73
74 TV_0 = 1; // D-wave K*2(1430)
75 T1_0 = 1;
76 T2_0 = 1;
77 m0_1430 = 1.4324;
78 width0_1430 = 0.109;
79 rho_1430 = 15;
80 phi_1430 = 0;
81 type[3] = 3;
82
83 mD = 1.96834; //Ds Mass
84 mPi = 0.49368; //K+ Mass
85 mK = 0.49368; //K- Mass
86
87 Pi = atan2(0.0,-1.0);
88 root2 = sqrt(2.);
89 root2d3 = sqrt(2./3);
90 root1d2 = sqrt(0.5);
91 root3d2 = sqrt(1.5);
92}
93
94void EvtDsToKKmunu::initProbMax() {
95 setProbMax(462526.1);
96}
97
98void EvtDsToKKmunu::decay( EvtParticle *p){
99/*
100 double maxprob = 0.0;
101 for(int ir=0;ir<=2000000;ir++){
102 p->initializePhaseSpace(getNDaug(),getDaugs());
103 EvtVector4R _K = p->getDaug(0)->getP4();
104 EvtVector4R _pi = p->getDaug(1)->getP4();
105 EvtVector4R _e = p->getDaug(2)->getP4();
106 EvtVector4R _nu = p->getDaug(3)->getP4();
107 int pid = EvtPDL::getStdHep(p->getDaug(0)->getId());
108 int charm;
109 if(pid == -321) charm = 1;
110 else charm = -1;
111 double m2, q2, cosV, cosL, chi;
112 KinVGen(_K, _pi, _e, _nu, charm, m2, q2, cosV, cosL, chi);
113 double _prob = calPDF(m2, q2, cosV, cosL, chi);
114 if(_prob>maxprob) {
115 maxprob=_prob;
116 std::cout << "Max PDF = " << ir << " charm= " << charm << " prob= " << _prob << std::endl;
117 }
118 }
119 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
120*/
121 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
122 ///////////////////////////////////////////////////////////////////////////////////////// test test test PDF ////////////////////////////////////////////////////////////////////////////////////////////////
123
124/*
125 p->initializePhaseSpace(getNDaug(),getDaugs());
126 EvtVector4R _K = p->getDaug(0)->getP4();
127 EvtVector4R _pi = p->getDaug(1)->getP4();
128 EvtVector4R _e = p->getDaug(2)->getP4();
129 EvtVector4R _nu = p->getDaug(3)->getP4();
130 int pid = EvtPDL::getStdHep(p->getDaug(0)->getId());
131 int charm;
132 if(pid == -321) charm = 1;
133 else charm = -1;
134
135 double m2, q2, cosV, cosL, chi;
136 KinVGen(_K, _pi, _e, _nu, charm, m2, q2, cosV, cosL, chi);
137 //m2 = 0.7317376201248132;
138 //q2 = 0.7930526827833025;
139 //cosV = -0.1044604580373859;
140 //cosL = 0.0503218879658880;
141 //chi = 1.9753715152960665;
142 //double _prob = calPDF(m2, q2, cosV, cosL, chi);
143
144 double _prob = calPDF(0.7317376201248132, 0.7930526827833025, -0.1044604580373859, 0.0503218879658880, 1.9753715152960665);
145 std::cout << "PDF PDF1 !!!!!!!!!! = " << _prob << std::endl;
146 _prob = calPDF(0.8253460062251410, 0.7590531075741469, -0.5990247196671010, -0.9155344197677775, -1.3296895383709690);
147 std::cout << "PDF PDF2 !!!!!!!!!! = " << _prob << std::endl;
148 _prob = calPDF(0.8161884101984322, 0.2477364661970792, -0.2764898614626146, 0.2826827644931852, 0.9307086619444224);
149 std::cout << "PDF PDF3 !!!!!!!!!! = " << _prob << std::endl;
150 _prob = calPDF(1.0912051424949705, 0.1906228414935633, 0.8132818897807017, -0.1032400648183377, 0.9686208243587908);
151 std::cout << "PDF PDF4 !!!!!!!!!! = " << _prob << std::endl;
152 _prob = calPDF(0.8079301573422222, 0.9019596037894472, 0.7789653251154771, 0.5324409155087387, -2.4926364105516909);
153 std::cout << "PDF PDF5 !!!!!!!!!! = " << _prob << std::endl;
154 _prob = calPDF(0.6698684771433824, 0.1744218225194828, 0.7286179682544673, -0.1705130677785587, 1.8691658148034282);
155 std::cout << "PDF PDF6 !!!!!!!!!! = " << _prob << std::endl;
156
157 double _prob = calPDF(1.0754665082688288, 0.6798367607443577, -0.1726785245179900, -0.8252479196790240, 2.8955919246019519);
158 std::cout << "PDF PDF1 !!!!!!!!!! = " << _prob << std::endl;
159 _prob = calPDF(1.0386029533289431, 0.4389802155039476, -0.5990247196671010, -0.9155344197677775, -1.3296895383709690);
160 std::cout << "PDF PDF2 !!!!!!!!!! = " << _prob << std::endl;
161 _prob = calPDF(1.0307383656072204, 0.2894701649730593, -0.2764898614626146, 0.2826827644931852, 0.9307086619444224);
162 std::cout << "PDF PDF3 !!!!!!!!!! = " << _prob << std::endl;
163 _prob = calPDF(1.0501028363282328, 0.1987927934599780, 0.8132818897807017, -0.1032400648183377, 0.9686208243587908);
164 std::cout << "PDF PDF4 !!!!!!!!!! = " << _prob << std::endl;
165 _prob = calPDF(1.0433160107621029, 0.2206250677498859, 0.7789653251154771, 0.5324409155087387, -2.4926364105516909);
166 std::cout << "PDF PDF5 !!!!!!!!!! = " << _prob << std::endl;
167 _prob = calPDF(1.0493574400231911, 0.2460442691281926, 0.7286179682544673, -0.1705130677785587, 1.8691658148034282);
168
169*/
170 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
171 /////////////////////////////////////////////////////////////////////////////////////////////// the end //////////////////////////////////////////////////////////////////////////////////////////////////////
172
173 p->initializePhaseSpace(getNDaug(),getDaugs());
174 EvtVector4R K = p->getDaug(0)->getP4();
175 EvtVector4R pi = p->getDaug(1)->getP4();
176 EvtVector4R e = p->getDaug(2)->getP4();
177 EvtVector4R nu = p->getDaug(3)->getP4();
178
179 int pid = EvtPDL::getStdHep(p->getDaug(0)->getId());
180 int charm;
181 if(pid == -321) charm = 1;
182 else charm = -1;
183 double m2, q2, cosV, cosL, chi;
184 KinVGen(K, pi, e, nu, charm, m2, q2, cosV, cosL, chi);
185 double prob = calPDF(m2, q2, cosV, cosL, chi);
186 setProb(prob);
187 return;
188
189}
190
191void EvtDsToKKmunu::KinVGen(EvtVector4R vp4_K, EvtVector4R vp4_Pi, EvtVector4R vp4_Lep, EvtVector4R vp4_Nu, int charm, double& m2, double& q2, double& cosV, double& cosL, double& chi)
192{
193 EvtVector4R vp4_KPi = vp4_K + vp4_Pi;
194 EvtVector4R vp4_W = vp4_Lep + vp4_Nu;
195
196 m2 = vp4_KPi.mass2();
197 q2 = vp4_W.mass2();
198
199 EvtVector4R boost;
200 boost.set(vp4_KPi.get(0), -vp4_KPi.get(1), -vp4_KPi.get(2), -vp4_KPi.get(3));
201 EvtVector4R vp4_Kp = boostTo(vp4_K, boost);
202 cosV = vp4_Kp.dot(vp4_KPi)/(vp4_Kp.d3mag()*vp4_KPi.d3mag());
203
204 boost.set(vp4_W.get(0), -vp4_W.get(1), -vp4_W.get(2), -vp4_W.get(3));
205 EvtVector4R vp4_Lepp = boostTo(vp4_Lep, boost);
206 cosL = vp4_Lepp.dot(vp4_W)/(vp4_Lepp.d3mag()*vp4_W.d3mag());
207
208 EvtVector4R V = vp4_KPi/vp4_KPi.d3mag();
209 EvtVector4R C = vp4_Kp.cross(V);
210 C/=C.d3mag();
211 EvtVector4R D = vp4_Lepp.cross(V);
212 D/=D.d3mag();
213 double sinx = C.cross(V).dot(D);
214 double cosx = C.dot(D);
215 chi = sinx>0? acos(cosx): -acos(cosx);
216 if(charm==-1) chi=-chi;
217}
218
219double EvtDsToKKmunu::calPDF(double m2, double q2, double cosV, double cosL, double chi) {
220 double delta[5] = {0};
221 double amplitude[5] = {0};
222 double m = sqrt(m2);
223 double q = sqrt(q2);
224
225 //begin to calculate form factor
226 EvtComplex F10(0.0,0.0);
227 EvtComplex F20(0.0,0.0);
228 EvtComplex F30(0.0,0.0);
229 EvtComplex F40(0.0,0.0);
230 EvtComplex F11(0.0,0.0);
231 EvtComplex F21(0.0,0.0);
232 EvtComplex F31(0.0,0.0);
233 EvtComplex F41(0.0,0.0);
234 EvtComplex F12(0.0,0.0);
235 EvtComplex F22(0.0,0.0);
236 EvtComplex F32(0.0,0.0);
237 EvtComplex F42(0.0,0.0);
238
239 EvtComplex f10(0.0,0.0);
240 EvtComplex f20(0.0,0.0);
241 EvtComplex f30(0.0,0.0);
242 EvtComplex f40(0.0,0.0);
243 EvtComplex f11(0.0,0.0);
244 EvtComplex f21(0.0,0.0);
245 EvtComplex f31(0.0,0.0);
246 EvtComplex f41(0.0,0.0);
247 EvtComplex f12(0.0,0.0);
248 EvtComplex f22(0.0,0.0);
249 EvtComplex f32(0.0,0.0);
250 EvtComplex f42(0.0,0.0);
251 EvtComplex coef(0.0,0.0);
252 double amplitude_temp, delta_temp;
253
254 for(int index=1; index<nAmps; index++)
255 {
256 switch(type[index])
257 {
258 case 0: //calculate form factor of S wave
259 {
260 NRS(m, q, rS, rS1, a_delta, b_delta, mA, m0_1430_S, width0_1430_S, amplitude_temp, delta_temp, f10,f40);
261 amplitude[index] = amplitude_temp;
262 delta[index] = delta_temp;
263 F10 = F10 + f10;
264 F40 = F40 + f40;
265 break;
266 }
267 case 1: //calculate form factor of P wave (phi)
268 {
269 ResonanceP(m, q, mV, mA, V_0, A1_0, A2_0, m0, width0, rBW, amplitude_temp, delta_temp, f11, f21, f31,f41);
270 amplitude[index] = amplitude_temp;
271 delta[index] = delta_temp;
272 coef = getCoef(rho, phi);
273 F11 = F11+ coef*f11;
274 F21 = F21+ coef*f21;
275 F31 = F31+ coef*f31;
276 F41 = F41+ coef*f41;
277 break;
278 }
279 case 2: //calculate form factor of P wave (K*(1410))
280 {
281
282 ResonanceP(m, q, mV, mA, V_0, A1_0, A2_0, m0_1410, width0_1410, rBW, amplitude_temp, delta_temp, f11, f21, f31,f41);
283 amplitude[index] = amplitude_temp;
284 delta[index] = delta_temp;
285 coef = getCoef(rho_1410, phi_1410);
286 F11 = F11+ coef*f11;
287 F21 = F21+ coef*f21;
288 F31 = F31+ coef*f31;
289 F41 = F41+ coef*f41;
290 break;
291 }
292 case 3: //calculate form factor of D wave
293 {
294 ResonanceD(m, q, mV, mA, TV_0, T1_0, T2_0, m0_1430, width0_1430, rBW, amplitude_temp, delta_temp, f12, f22, f32);
295 amplitude[index] = amplitude_temp;
296 delta[index] = delta_temp;
297 coef = getCoef(rho_1430, phi_1430);
298 F12 = F12+ coef*f12;
299 F22 = F22+ coef*f22;
300 F32 = F32+ coef*f32;
301 break;
302 }
303 default:
304 {
305 std::cout<<"No such form factor type!!!"<<std::endl;
306 break;
307 }
308 }
309 }
310
311 //begin to calculate pdf value
312 double I,I1,I2,I3,I4,I5,I6,I7,I8,I9;
313
314 double cosV2 = cosV*cosV;
315 double sinV = sqrt(1.0-cosV2);
316 double sinV2 = sinV*sinV;
317 double mu=0.105658;
318 double betaL =(1- mu*mu/q2);
319
320 EvtComplex F1 = F10 + F11*cosV + F12*(1.5*cosV2-0.5);
321 EvtComplex F2 = F21*root1d2 + F22*cosV*root3d2;
322 EvtComplex F3 = F31*root1d2 + F32*cosV*root3d2;
323 EvtComplex F4 = F40+F41*(cosV);
324
325// I1 = 0.25*( abs2(F1) + 1.5*sinV2*(abs2(F2) + abs2(F3) ) );
326// I2 = -0.25*( abs2(F1) - 0.5*sinV2*(abs2(F2) + abs2(F3) ) );
327// I3 = -0.25*( abs2(F2) - abs2(F3) )*sinV2;
328// I4 = real( conj(F1)*F2 )*sinV*0.5;
329// I5 = real( conj(F1)*F3 )*sinV;
330// I6 = real( conj(F2)*F3 )*sinV2;
331// I7 = imag( conj(F2)*F1 )*sinV;
332// I8 = imag( conj(F3)*F1 )*sinV*0.5;
333// I9 = imag( conj(F3)*F2 )*sinV2*(-0.5);
334
335 I1 = 0.25*(2- betaL)*betaL*real((F1*(conj(F1)))) + ((betaL/2.0)-(betaL*betaL/8.0))*(sinV2)*(real((F2*conj(F2))) + real((F3*conj(F3)))) + 0.5*(1- betaL)*betaL*real((F4*(conj(F4))));
336 I2 = -1/4.0*( betaL*betaL)*(real((F1*conj(F1)))-0.5*sinV2*(real((F2*conj(F2))) + real((F3*conj(F3)))));
337 I3 = -0.25*(betaL*betaL)*(real((F2*conj(F2))) - real((F3*conj(F3))))*sinV2;
338 I4 = real((conj(F2)*F1))*sinV*0.5*(betaL*betaL);
339 I5 = real((conj(F3)*F1 - ((1-betaL))*conj(F4)*F2))*sinV*(betaL);
340 I6 = (betaL)*real((conj(F3)*F2*sinV2 + (1-betaL)*conj(F4)*F1));
341 I7 = imag((conj(F2)*F1*sinV*(betaL)+sinV*(betaL)*(1-betaL)*conj(F4)*F3));
342 I8 = imag((conj(F3)*F1))*sinV*(0.5)*(betaL*betaL);
343 I9 = imag((conj(F3)*F2))*sinV2*(-0.5)*(betaL*betaL);
344
345// I1 = 0.25*(2- betaL)*betaL*(F1*(std::conj(F1))).real() + ((betaL/2.0)-(betaL*betaL/8.0))*(sinV2)*((F2*std::conj(F2)).real() + (F3*std::conj(F3)).real()) + 0.5*(1- betaL)*betaL*(F4*(std::conj(F4))).real();
346 // I2 = -1/4.0*( betaL*betaL)*((F1*std::conj(F1)).real() -0.5*sinV2*((F2*std::conj(F2)).real() + (F3*std::conj(F3)).real()));
347 // I3 = -0.25*(betaL*betaL)*((F2*std::conj(F2)).real() - (F3*std::conj(F3)).real())*sinV2;
348 // I4 = (std::conj(F2)*F1).real()*sinV*0.5*(betaL*betaL);
349 // I5 = (std::conj(F3)*F1 - ((1-betaL))*std::conj(F4)*F2).real()*sinV*(betaL);
350 // I6 = (betaL)*(std::conj(F3)*F2*sinV2 + (1-betaL)*std::conj(F4)*F1).real();
351 // I7 = (std::conj(F2)*F1*sinV*(betaL)+sinV*(betaL)*(1-betaL)*std::conj(F4)*F3).imag();
352 // I8 = (std::conj(F3)*F1).imag()*sinV*(0.5)*(betaL*betaL);
353 // I9 = (std::conj(F3)*F2).imag()*sinV2*(-0.5)*(betaL*betaL);
354
355
356
357 double coschi = cos(chi);
358 double sinchi = sin(chi);
359 double sin2chi = 2.0*sinchi*coschi;
360 double cos2chi = 1.0 - 2.0*sinchi*sinchi;
361
362 double sinL = sqrt(1.-cosL*cosL);
363 double sinL2 = sinL*sinL;
364 double sin2L = 2.0*sinL*cosL;
365 double cos2L = 1.0 - 2.0*sinL2;
366
367 I = I1 + I2*cos2L + I3*sinL2*cos2chi + I4*sin2L*coschi + I5*sinL*coschi + I6*cosL + I7*sinL*sinchi + I8*sin2L*sinchi + I9*sinL2*sin2chi;
368 return I;
369}
370
371void EvtDsToKKmunu::ResonanceP(double m, double q, double mV, double mA, double V_0, double A1_0, double A2_0, double m0, double width0, double rBW, double& amplitude, double& delta, EvtComplex& F11, EvtComplex& F21, EvtComplex& F31, EvtComplex& F41)
372{
373 double pKPi = getPStar(mD, m, q);
374 double mD2 = mD*mD;
375 double m2 = m*m;
376 double m02 = m0*m0;
377 double q2 = q*q;
378 double mV2 = mV*mV;
379 double mA2 = mA*mA;
380 double summDm = mD+m;
381 double V = V_0 /(1.0-q2/(mV2));
382 double A1 = A1_0/(1.0-q2/(mA2));
383 double A2 = A2_0/(1.0-q2/(mA2));
384 double A = summDm*A1;
385 double B = 2.0*mD*pKPi/summDm*V;
386 double A0 = ((mD+m)*A1-(mD-m)*A2)/(2*m);
387
388 //construct the helicity form factor
389 double H0 = 0.5/(m*q)*((mD2-m2-q2)*summDm*A1-4.0*(mD2*pKPi*pKPi)/summDm*A2);
390 double Hp = A-B;
391 double Hm = A+B;
392
393 //calculate alpha
394 // double B_Kstar = 1./3.; //B_Kstar = Br(Kstar(892)->k pi)
395 double B_Kstar = 0.4920; //B_Phi = Br(Phi->k k)
396 double pStar0 = getPStar(m0, mPi, mK);
397 double alpha = sqrt(3.*Pi*B_Kstar/(pStar0*width0));
398
399 //construct amplitudes of (non)resonance
400 double F = getF1(m, m0, mPi, mK, rBW);
401 double width = getWidth1(m, m0, mPi, mK, width0, rBW);
402
403 EvtComplex C(m0*width0*F,0.0);
404 double AA = m02-m2;
405 double BB = -m0*width;
406 EvtComplex amp = C/EvtComplex(AA,BB);
407 amplitude = abs(amp);
408 delta = atan2(imag(amp), real(amp));
409
410 double alpham2 = alpha*2.0;
411 double alpham3 = alpha*4.0;
412 F11 = amp*alpham2*q*H0*root2;
413 F21 = amp*alpham2*q*(Hp+Hm);
414 F31 = amp*alpham2*q*(Hp-Hm);
415 F41 = amp*alpham3*q*(mD*pKPi*A0/q)*root2;
416}
417
418void EvtDsToKKmunu::NRS(double m, double q, double rS, double rS1, double a_delta, double b_delta, double mA, double m0, double width0, double& amplitude, double& delta, EvtComplex& F10,EvtComplex& F40)
419{
420 static const double tmp = (mK+mPi)*(mK+mPi);
421
422 double m2 = m*m;
423 double q2 = q*q;
424 double mA2 = mA*mA;
425 double pKPi = getPStar(mD, m, q);
426 double m_K0_1430 = m0;
427 double width_K0_1430 = width0;
428 double m2_K0_1430 = m_K0_1430*m_K0_1430;
429 double width = getWidth0(m, m_K0_1430, mPi, mK, width_K0_1430);
430
431 //calculate modul of the amplitude
432 double x,Pm;
433 if(m<m_K0_1430) {
434 x = sqrt(m2/tmp-1.0);
435 Pm = 1.0+rS1*x;
436 } else {
437 x = sqrt(m2_K0_1430/tmp-1.0);
438 Pm = 1.0+rS1*x;
439 Pm *= m_K0_1430*width_K0_1430/sqrt((m2_K0_1430-m2)*(m2_K0_1430-m2)+m2_K0_1430*width*width);
440 }
441
442 //calculate phase of the amplitude
443 double pStar = getPStar(m, mPi, mK);
444 double delta_bg = atan(2.*a_delta*pStar/(2.+a_delta*b_delta*pStar*pStar));
445 delta_bg = (delta_bg>0)? delta_bg: (delta_bg+Pi);
446
447 double delta_K0_1430 = atan(m_K0_1430*width/(m2_K0_1430-m2));
448 delta_K0_1430 = (delta_K0_1430>0)? delta_K0_1430: (delta_K0_1430+Pi);
449 delta = delta_bg + delta_K0_1430;
450
451 EvtComplex ci(cos(delta), sin(delta));
452 EvtComplex amp = ci*rS*Pm;
453 amplitude = rS*Pm;
454
455 F10 = amp*pKPi*mD/(1.-q2/mA2);
456 F40 = amp*q2/(1-q2/mA2);
457}
458
459void EvtDsToKKmunu::ResonanceD(double m, double q, double mV, double mA, double TV_0, double T1_0, double T2_0, double m0, double width0, double rBW, double& amplitude, double& delta, EvtComplex& F12, EvtComplex& F22, EvtComplex& F32)
460{
461 double pKPi = getPStar(mD, m, q);
462 double mD2 = mD*mD;
463 double m2 = m*m;
464 double m02 = m0*m0;
465 double q2 = q*q;
466 double mV2 = mV*mV;
467 double mA2 = mA*mA;
468 double summDm = mD+m;
469 double TV = TV_0 /(1.0-q2/(mV2));
470 double T1 = T1_0/(1.0-q2/(mA2));
471 double T2 = T2_0/(1.0-q2/(mA2));
472
473 //construct amplitudes of (non)resonance
474 double F = getF2(m, m0, mPi, mK, rBW);
475 double width = getWidth2(m, m0, mPi, mK, width0, rBW);
476 EvtComplex C(m0*width0*F,0.0);
477 double AA = m02-m2;
478 double BB = -m0*width;
479
480 EvtComplex amp = C/EvtComplex(AA,BB);
481 amplitude = abs(amp);
482 delta = atan2(imag(amp), real(amp));
483
484 F12 = amp*mD*pKPi/3.*((mD2-m2-q2)*summDm*T1-mD2*pKPi*pKPi/summDm*T2);
485 F22 = amp*root2d3*mD*m*q*pKPi*summDm*T1;
486 F32 = amp*root2d3*2.*mD2*m*q*pKPi*pKPi/summDm*TV;
487}
488
489double EvtDsToKKmunu::getPStar(double m, double m1, double m2)
490{
491 double s = m*m;
492 double s1 = m1*m1;
493 double s2 = m2*m2;
494 double x = s + s1 - s2;
495 double t = 0.25*x*x/s - s1;
496 double p;
497 if (t>0.0) {
498 p = sqrt(t);
499 } else {
500 std::cout << " Hello, pstar is less than 0.0" << std::endl;
501 p = 0.04;
502 }
503 return p;
504}
505
506double EvtDsToKKmunu::getF1(double m, double m0, double m_c1, double m_c2, double rBW)
507{
508 double pStar = getPStar(m, m_c1, m_c2);
509 double pStar0 = getPStar(m0, m_c1, m_c2);
510 double rBW2 = rBW*rBW;
511 double pStar2 = pStar*pStar;
512 double pStar02 = pStar0*pStar0;
513 double B = 1./sqrt(1.+rBW2*pStar2);
514 double B0 = 1./sqrt(1.+rBW2*pStar02);
515 double F = pStar/pStar0*B/B0;
516 return F;
517}
518
519double EvtDsToKKmunu::getF2(double m, double m0, double m_c1, double m_c2, double rBW)
520{
521 double pStar = getPStar(m, m_c1, m_c2);
522 double pStar0 = getPStar(m0, m_c1, m_c2);
523 double rBW2 = rBW*rBW;
524 double pStar2 = pStar*pStar;
525 double pStar02 = pStar0*pStar0;
526 double B = 1./sqrt((rBW2*pStar2-3.)*(rBW2*pStar2-3.)+9.*rBW2*pStar2);
527 double B0 = 1./sqrt((rBW2*pStar02-3.)*(rBW2*pStar02-3.)+9.*rBW2*pStar02);
528 double F = pStar2/pStar02*B/B0;
529 return F;
530}
531
532double EvtDsToKKmunu::getWidth0(double m, double m0, double m_c1, double m_c2, double width0)
533{
534 double pStar = getPStar(m, m_c1, m_c2);
535 double pStar0 = getPStar(m0, m_c1, m_c2);
536 double width = width0*pStar/pStar0*m0/m;
537 return width;
538}
539
540double EvtDsToKKmunu::getWidth1(double m, double m0, double m_c1, double m_c2, double width0, double rBW)
541{
542 double pStar = getPStar(m, m_c1, m_c2);
543 double pStar0 = getPStar(m0, m_c1, m_c2);
544 double F = getF1(m, m0, m_c1, m_c2, rBW);
545 double width = width0*pStar/pStar0*m0/m*F*F;
546 return width;
547}
548
549double EvtDsToKKmunu::getWidth2(double m, double m0, double m_c1, double m_c2, double width0, double rBW)
550{
551 double pStar = getPStar(m, m_c1, m_c2);
552 double pStar0 = getPStar(m0, m_c1, m_c2);
553 double F = getF2(m, m0, m_c1, m_c2, rBW);
554 double width = width0*pStar/pStar0*m0/m*F*F;
555 return width;
556}
557
558EvtComplex EvtDsToKKmunu::getCoef(double rho, double phi)
559{
560 EvtComplex coef ( rho*cos(phi), rho*sin(phi) );
561 return coef;
562}
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
delta
const double delta
Definition CheckRung/check/gainVsDoca.cxx:32
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
I
const DifComplex I
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition Evt3Rank3C.cc:175
imag
double imag(const EvtComplex &c)
Definition EvtComplex.hh:246
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition EvtDiracSpinor.cc:75
alpha
const double alpha
Definition FastVertexFit.cxx:4
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
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29
t
TTree * t
Definition binning.cxx:23
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:533
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:520
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:297
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:504
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:65
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:482
EvtDecayProb::setProb
void setProb(double prob)
Definition EvtDecayProb.hh:34
EvtDsToKKmunu::decay
void decay(EvtParticle *p)
Definition EvtDsToKKmunu.cc:98
EvtDsToKKmunu::getName
void getName(std::string &name)
Definition EvtDsToKKmunu.cc:29
EvtDsToKKmunu::~EvtDsToKKmunu
virtual ~EvtDsToKKmunu()
Definition EvtDsToKKmunu.cc:27
EvtDsToKKmunu::clone
EvtDecayBase * clone()
Definition EvtDsToKKmunu.cc:33
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
EvtVector4R::dot
double dot(const EvtVector4R &v2) const
Definition EvtVector4R.cc:199
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:179
EvtVector4R::cross
EvtVector4R cross(const EvtVector4R &v2)
Definition EvtVector4R.cc:171
EvtVector4R::d3mag
double d3mag() const
Definition EvtVector4R.cc:186
EvtVector4R::mass2
double mass2() const
Definition EvtVector4R.hh:116
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:183
m1
double * m1
Definition qcdloop1.h:75
m2
double double * m2
Definition qcdloop1.h:75
pi
const float pi
Definition vector3.h:133