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