BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDTopipienu.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: EvtDTopipienu.cc
11// the necessary file: EvtDTopipienu.hh
12//
13// Description: D -> pi pi e+ nu,
14// PHYSICAL REVIEW Letter 122, 062001 (2019)
15//
16// Modification history:
17//
18// Liaoyuan Dong Jan.16, 2020 Module created
19//
20//------------------------------------------------------------------------
21#include "EvtDTopipienu.hh"
22#include "../EvtGenBase/EvtComplex.hh"
23#include "../EvtGenBase/EvtDecayTable.hh"
24#include "../EvtGenBase/EvtGenKine.hh"
25#include "../EvtGenBase/EvtPDL.hh"
26#include "../EvtGenBase/EvtParticle.hh"
27#include "../EvtGenBase/EvtPatches.hh"
28#include "../EvtGenBase/EvtReport.hh"
29#include <stdlib.h>
30
31EvtDTopipienu::~EvtDTopipienu() {}
32
33void EvtDTopipienu::getName( std::string& model_name ) { model_name = "DTopipienu"; }
34
35EvtDecayBase* EvtDTopipienu::clone() { return new EvtDTopipienu; }
36
37void EvtDTopipienu::init() {
38 checkNArg( 0 );
39 checkNDaug( 4 );
40 checkSpinParent( EvtSpinType::SCALAR );
41 checkSpinDaughter( 0, EvtSpinType::SCALAR );
42 checkSpinDaughter( 1, EvtSpinType::SCALAR );
43
44 first = 1;
45 last = 4;
46 ProbMax = 591000;
47 // cout << "Initializing EvtDTopipienu: ProbMax = " << ProbMax << endl;
48
49 type[0] = 0;
50 type[1] = 1;
51 type[2] = 2;
52 type[3] = 3;
53
54 mV = 2.01;
55 mA = 2.42;
56 V_0 = 1.6948;
57 A1_0 = 1;
58 A2_0 = 0.84489;
59
60 m0 = 0.77526;
61 width0 = 0.14910;
62 rBW = 3.0;
63 rho = 1.0;
64 phi = 0.0;
65 BF = 1.0;
66
67 m0_omega = 0.78265;
68 width0_omega = 0.00849;
69 rho_omega = 0.12902;
70 phi_omega = 2.9285;
71 BF_omega = 1.0;
72
73 m0_S = 0.953;
74 rho_S = 135.27;
75 phi_S = 3.4044;
76
77 Dp_mD = 1.86962;
78 Dp_mPi1 = 0.13957;
79 Dp_mPi2 = 0.13957;
80 D0_mD = 1.86486;
81 D0_mPi1 = 0.13957;
82 D0_mPi2 = 0.1349766;
83
84 Pi = atan2( 0.0, -1.0 );
85 root2 = sqrt( 2. );
86 root2d3 = sqrt( 2. / 3 );
87 root1d2 = sqrt( 0.5 );
88 root3d2 = sqrt( 1.5 );
89
90 mKa = 0.493677;
91 mPi = 0.13957;
92 mEt = 0.547853;
93}
94
95void EvtDTopipienu::initProbMax() { setProbMax( ProbMax ); }
96
97void EvtDTopipienu::decay( EvtParticle* p ) {
98 /*
99 double maxprob = 0.0;
100 for(int ir=0;ir<=60000000;ir++){
101 p->initializePhaseSpace(getNDaug(),getDaugs());
102 EvtVector4R _pi1 = p->getDaug(0)->getP4();
103 EvtVector4R _pi2 = p->getDaug(1)->getP4();
104 EvtVector4R _e = p->getDaug(2)->getP4();
105 EvtVector4R _nu = p->getDaug(3)->getP4();
106
107 int pid = EvtPDL::getStdHep(p->getDaug(0)->getId());
108 int charm;
109 if(pid == -211) charm = 1;
110 else charm = -1;
111 double m2, q2, cosV, cosL, chi;
112 KinVGen(_pi1, _pi2, _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 <<
117 std::endl;
118 }
119 }
120 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
121 */
122 p->initializePhaseSpace( getNDaug(), getDaugs() );
123 EvtVector4R pi1 = p->getDaug( 0 )->getP4();
124 EvtVector4R pi2 = p->getDaug( 1 )->getP4();
125 EvtVector4R e = p->getDaug( 2 )->getP4();
126 EvtVector4R nu = p->getDaug( 3 )->getP4();
127
128 int pid = EvtPDL::getStdHep( p->getDaug( 0 )->getId() );
129 int charm;
130 if ( pid == -211 ) charm = 1;
131 else charm = -1;
132 double m2, q2, cosV, cosL, chi;
133 KinVGen( pi1, pi2, e, nu, charm, m2, q2, cosV, cosL, chi );
134 double prob = calPDF( m2, q2, cosV, cosL, chi );
135 setProb( prob );
136 return;
137}
138
139void EvtDTopipienu::KinVGen( EvtVector4R vp4_K, EvtVector4R vp4_Pi, EvtVector4R vp4_Lep,
140 EvtVector4R vp4_Nu, int charm, double& m2, double& q2,
141 double& cosV, double& cosL, double& chi ) {
142 EvtVector4R vp4_KPi = vp4_K + vp4_Pi;
143 EvtVector4R vp4_W = vp4_Lep + vp4_Nu;
144 m2 = vp4_KPi.mass2();
145 q2 = vp4_W.mass2();
146
147 EvtVector4R boost;
148 boost.set( vp4_KPi.get( 0 ), -vp4_KPi.get( 1 ), -vp4_KPi.get( 2 ), -vp4_KPi.get( 3 ) );
149 EvtVector4R vp4_Kp = boostTo( vp4_K, boost );
150 cosV = vp4_Kp.dot( vp4_KPi ) / ( vp4_Kp.d3mag() * vp4_KPi.d3mag() );
151
152 boost.set( vp4_W.get( 0 ), -vp4_W.get( 1 ), -vp4_W.get( 2 ), -vp4_W.get( 3 ) );
153 EvtVector4R vp4_Lepp = boostTo( vp4_Lep, boost );
154 cosL = vp4_Lepp.dot( vp4_W ) / ( vp4_Lepp.d3mag() * vp4_W.d3mag() );
155
156 EvtVector4R V = vp4_KPi / vp4_KPi.d3mag();
157 EvtVector4R C = vp4_Kp.cross( V );
158 C /= C.d3mag();
159 EvtVector4R D = vp4_Lepp.cross( V );
160 D /= D.d3mag();
161 double sinx = C.cross( V ).dot( D );
162 double cosx = C.dot( D );
163 chi = sinx > 0 ? acos( cosx ) : -acos( cosx );
164 if ( charm == -1 ) chi = -chi;
165}
166
167double EvtDTopipienu::calPDF( double m2, double q2, double cosV, double cosL, double chi ) {
168 double m = sqrt( m2 );
169 double q = sqrt( q2 );
170
171 // begin to calculate form factor
172 EvtComplex F10( 0.0, 0.0 );
173 EvtComplex F11( 0.0, 0.0 );
174 EvtComplex F21( 0.0, 0.0 );
175 EvtComplex F31( 0.0, 0.0 );
176 EvtComplex F12( 0.0, 0.0 );
177 EvtComplex F22( 0.0, 0.0 );
178 EvtComplex F32( 0.0, 0.0 );
179
180 EvtComplex f10( 0.0, 0.0 );
181 EvtComplex f11( 0.0, 0.0 );
182 EvtComplex f21( 0.0, 0.0 );
183 EvtComplex f31( 0.0, 0.0 );
184
185 EvtComplex coef( 0.0, 0.0 );
186
187 for ( int index = first; index < last; index++ )
188 {
189 switch ( type[index] )
190 {
191 case 0: // calculate form factor of P wave (rho)
192 {
193 ResonanceGS( m, q, D0_mD, D0_mPi1, D0_mPi2, f11, f21, f31 );
194 coef = getCoef( rho, phi );
195 F11 = F11 + coef * f11;
196 F21 = F21 + coef * f21;
197 F31 = F31 + coef * f31;
198 break;
199 }
200 case 1: // calculate form factor of P wave (GS)
201 {
202 ResonanceGS( m, q, Dp_mD, Dp_mPi1, Dp_mPi2, f11, f21, f31 );
203 coef = getCoef( rho, phi );
204 F11 = F11 + coef * f11;
205 F21 = F21 + coef * f21;
206 F31 = F31 + coef * f31;
207 break;
208 }
209 case 2: // calculate form factor of omega
210 {
211 ResonancePGScbw( m, q, f11, f21, f31 );
212 coef = getCoef( rho_omega, phi_omega );
213 F11 = F11 + coef * f11;
214 F21 = F21 + coef * f21;
215 F31 = F31 + coef * f31;
216 break;
217 }
218 case 3: // calculate form factor of S wave (BW corrected by Bugg)
219 {
220 ResonanceSBugg( m, q, f10 );
221 coef = getCoef( rho_S, phi_S );
222 F10 = F10 + coef * f10;
223 break;
224 }
225 default: {
226 std::cout << "No such form factor type!!!" << std::endl;
227 break;
228 }
229 }
230 }
231
232 // begin to calculate pdf value
233 double I, I1, I2, I3, I4, I5, I6, I7, I8, I9;
234
235 double cosV2 = cosV * cosV;
236 double sinV = sqrt( 1.0 - cosV2 );
237 double sinV2 = sinV * sinV;
238
239 EvtComplex F1 = F10 + F11 * cosV + F12 * ( 1.5 * cosV2 - 0.5 );
240 EvtComplex F2 = F21 * root1d2 + F22 * cosV * root3d2;
241 EvtComplex F3 = F31 * root1d2 + F32 * cosV * root3d2;
242
243 I1 = 0.25 * ( abs2( F1 ) + 1.5 * sinV2 * ( abs2( F2 ) + abs2( F3 ) ) );
244 I2 = -0.25 * ( abs2( F1 ) - 0.5 * sinV2 * ( abs2( F2 ) + abs2( F3 ) ) );
245 I3 = -0.25 * ( abs2( F2 ) - abs2( F3 ) ) * sinV2;
246 I4 = real( conj( F1 ) * F2 ) * sinV * 0.5;
247 I5 = real( conj( F1 ) * F3 ) * sinV;
248 I6 = real( conj( F2 ) * F3 ) * sinV2;
249 I7 = imag( conj( F2 ) * F1 ) * sinV;
250 I8 = imag( conj( F3 ) * F1 ) * sinV * 0.5;
251 I9 = imag( conj( F3 ) * F2 ) * sinV2 * ( -0.5 );
252
253 double coschi = cos( chi );
254 double sinchi = sin( chi );
255 double sin2chi = 2.0 * sinchi * coschi;
256 double cos2chi = 1.0 - 2.0 * sinchi * sinchi;
257
258 double sinL = sqrt( 1. - cosL * cosL );
259 double sinL2 = sinL * sinL;
260 double sin2L = 2.0 * sinL * cosL;
261 double cos2L = 1.0 - 2.0 * sinL2;
262
263 I = I1 + I2 * cos2L + I3 * sinL2 * cos2chi + I4 * sin2L * coschi + I5 * sinL * coschi +
264 I6 * cosL + I7 * sinL * sinchi + I8 * sin2L * sinchi + I9 * sinL2 * sin2chi;
265 return I;
266}
267
268void EvtDTopipienu::ResonanceGS( double m, double q, double massD, double massPi1,
269 double massPi2, EvtComplex& F11, EvtComplex& F21,
270 EvtComplex& F31 ) {
271 double pKPi = getPStar( massD, m, q );
272 double mD2 = massD * massD;
273 double m2 = m * m;
274 double m02 = m0 * m0;
275 double q2 = q * q;
276 double mV2 = mV * mV;
277 double mA2 = mA * mA;
278 double summDm = massD + m;
279 double V = V_0 / ( 1.0 - q2 / ( mV2 ) );
280 double A1 = A1_0 / ( 1.0 - q2 / ( mA2 ) );
281 double A2 = A2_0 / ( 1.0 - q2 / ( mA2 ) );
282 double A = summDm * A1;
283 double B = 2.0 * massD * pKPi / summDm * V;
284
285 // construct the helicity form factor
286 double H0 = 0.5 / ( m * q ) *
287 ( ( mD2 - m2 - q2 ) * summDm * A1 - 4.0 * ( mD2 * pKPi * pKPi ) / summDm * A2 );
288 double Hp = A - B;
289 double Hm = A + B;
290
291 // calculate alpha
292 // double B_Kstar = 2./3.0; //B_Kstar = Br(Kstar(892)->k pi)
293 double pStar = getPStar( m, massPi1, massPi2 );
294 double pStar0 = getPStar( m0, massPi1, massPi2 );
295 double alpha = sqrt( 3.0 * Pi * BF / ( pStar0 * width0 ) );
296
297 // construct amplitudes of (non)resonance
298 double F = getF1( m, m0, massPi1, massPi2, rBW );
299 EvtComplex C( m02 * ( 1.0 + width0 * getGx( m0, pStar0, massPi1, massPi2 ) / m0 ) * F, 0.0 );
300 double AA = m02 - m2 + width0 * getFx( m02, m2, pStar, pStar0, massPi1, massPi2 );
301 double BB = -m0 * getWidthrho( m, m0, width0, pStar, pStar0 );
302 EvtComplex tmp( AA, BB );
303 EvtComplex amp = C / tmp;
304
305 double alpham2 = alpha * 2.0;
306 F11 = amp * alpham2 * q * H0 * root2;
307 F21 = amp * alpham2 * q * ( Hp + Hm );
308 F31 = amp * alpham2 * q * ( Hp - Hm );
309}
310
311void EvtDTopipienu::ResonancePGScbw( double m, double q, EvtComplex& F11, EvtComplex& F21,
312 EvtComplex& F31 ) {
313 double pKPi = getPStar( Dp_mD, m, q );
314 double mD2 = Dp_mD * Dp_mD;
315 double m2 = m * m;
316 double m02 = m0_omega * m0_omega;
317 double mR2 = m0 * m0;
318 double q2 = q * q;
319 double mV2 = mV * mV;
320 double mA2 = mA * mA;
321 double summDm = Dp_mD + m;
322 double V = V_0 / ( 1.0 - q2 / ( mV2 ) );
323 double A1 = A1_0 / ( 1.0 - q2 / ( mA2 ) );
324 double A2 = A2_0 / ( 1.0 - q2 / ( mA2 ) );
325 double A = summDm * A1;
326 double B = 2.0 * Dp_mD * pKPi / summDm * V;
327 // construct the helicity form factor
328 double H0 = 0.5 / ( m * q ) *
329 ( ( mD2 - m2 - q2 ) * summDm * A1 - 4.0 * ( mD2 * pKPi * pKPi ) / summDm * A2 );
330 double Hp = A - B;
331 double Hm = A + B;
332
333 // calculate alpha
334 // double B_Kstar = 2./3.0; //B_Kstar = Br(Kstar(892)->k pi)
335 double pStar = getPStar( m, Dp_mPi1, Dp_mPi2 );
336 double pStar0 = getPStar( m0_omega, Dp_mPi1, Dp_mPi2 );
337 double alpha = sqrt( 3.0 * Pi * BF_omega / ( pStar0 * width0_omega ) );
338
339 // construct amplitudes of (non)resonance
340 double F = getF1( m, m0_omega, Dp_mPi1, Dp_mPi2, rBW );
341 EvtComplex amp1( 0.0, 0.0 );
342 EvtComplex amp2( 0.0, 0.0 );
343 // CBW
344 EvtComplex C1( m0_omega * width0_omega * F, 0.0 );
345 double AA1 = m02 - m2;
346 double BB1 = -m0_omega * width0_omega;
347 EvtComplex tmp1( AA1, BB1 );
348 amp1 = C1 / tmp1;
349 // GS
350 pStar0 = getPStar( m0, Dp_mPi1, Dp_mPi2 );
351 EvtComplex C2( mR2 * ( 1.0 + width0 * getGx( m0, pStar0, Dp_mPi1, Dp_mPi2 ) / m0 ), 0.0 );
352 double AA2 = mR2 - m2 + width0 * getFx( mR2, m2, pStar, pStar0, Dp_mPi1, Dp_mPi2 );
353 double BB2 = -m0 * getWidthrho( m, m0, width0, pStar, pStar0 );
354 EvtComplex tmp2( AA2, BB2 );
355 amp2 = C2 / tmp2;
356 EvtComplex amp = amp1 * amp2;
357
358 double alpham2 = alpha * 2.0;
359 F11 = amp * alpham2 * q * H0 * root2;
360 F21 = amp * alpham2 * q * ( Hp + Hm );
361 F31 = amp * alpham2 * q * ( Hp - Hm );
362}
363
364void EvtDTopipienu::ResonanceSBugg( double m, double q, EvtComplex& F10 ) {
365 double pKPi = getPStar( Dp_mD, m, q );
366 double m2 = m * m;
367 double q2 = q * q;
368 double mA2 = mA * mA;
369
370 double sA = 0.41 * mPi * mPi;
371 double mr = m0_S; // 0.953;
372 double mr2 = m0_S * m0_S; // 0.908209;// 0.953*0.953;
373 double alpha = 1.3;
374 double g4pi = 0.011;
375
376 EvtComplex ciR( 1.0, 0.0 );
377 EvtComplex ciM( 0.0, 1.0 );
378 EvtComplex Gamma1( 0.0, 0.0 );
379 EvtComplex Gamma2( 0.0, 0.0 );
380 EvtComplex Gamma3( 0.0, 0.0 );
381 EvtComplex Gamma4( 0.0, 0.0 );
382
383 Gamma1 = getrho( m2, mPi ) * getG1( m2, mr ) * ( m2 - sA ) / ( mr2 - sA );
384 Gamma2 = getrho( m2, mKa ) * 0.6 * getG1( m2, mr ) * ( m2 / mr2 ) *
385 exp( -alpha * fabs( m2 - 4.0 * mKa * mKa ) );
386 Gamma3 = getrho( m2, mEt ) * 0.2 * getG1( m2, mr ) * ( m2 / mr2 ) *
387 exp( -alpha * fabs( m2 - 4.0 * mEt * mEt ) );
388 if ( m > 4 * mPi )
389 Gamma4 = ciR * mr * g4pi * ( 1.0 + exp( 7.082 - 2.845 * mr2 ) ) /
390 ( 1.0 + exp( 7.082 - 2.845 * m2 ) );
391
392 double AA = mr2 - m2 - getG1( m2, mr ) * ( m2 - sA ) * getZ( m2, mr2 ) / ( mr2 - sA );
393 EvtComplex amp = ciR / ( ciR * AA - ciM * ( Gamma1 + Gamma2 + Gamma3 + Gamma4 ) );
394 F10 = amp * pKPi * Dp_mD / ( 1.0 - q2 / mA2 );
395}
396
397double EvtDTopipienu::getPStar( double m, double m1, double m2 ) {
398 double s = m * m;
399 double s1 = m1 * m1;
400 double s2 = m2 * m2;
401 double x = s + s1 - s2;
402 double t = 0.25 * x * x / s - s1;
403 double p;
404 if ( t > 0.0 ) { p = sqrt( t ); }
405 else
406 {
407 // std::cout << " Hello, pstar is less than 0.0" << std::endl;
408 p = 0.04;
409 }
410 return p;
411}
412
413double EvtDTopipienu::getF1( double m, double m0, double m_c1, double m_c2, double rBW ) {
414 double pStar = getPStar( m, m_c1, m_c2 );
415 double pStar0 = getPStar( m0, m_c1, m_c2 );
416 double rBW2 = rBW * rBW;
417 double pStar2 = pStar * pStar;
418 double pStar02 = pStar0 * pStar0;
419 double B = 1. / sqrt( 1. + rBW2 * pStar2 );
420 double B0 = 1. / sqrt( 1. + rBW2 * pStar02 );
421 double F = pStar / pStar0 * B / B0;
422 return F;
423}
424
425double EvtDTopipienu::getF2( double m, double m0, double m_c1, double m_c2, double rBW ) {
426 double pStar = getPStar( m, m_c1, m_c2 );
427 double pStar0 = getPStar( m0, m_c1, m_c2 );
428 double rBW2 = rBW * rBW;
429 double pStar2 = pStar * pStar;
430 double pStar02 = pStar0 * pStar0;
431 double B = 1. / sqrt( ( rBW2 * pStar2 - 3. ) * ( rBW2 * pStar2 - 3. ) + 9. * rBW2 * pStar2 );
432 double B0 =
433 1. / sqrt( ( rBW2 * pStar02 - 3. ) * ( rBW2 * pStar02 - 3. ) + 9. * rBW2 * pStar02 );
434 double F = pStar2 / pStar02 * B / B0;
435 return F;
436}
437
438double EvtDTopipienu::getWidth0( double m, double m0, double m_c1, double m_c2,
439 double width0 ) {
440 double pStar = getPStar( m, m_c1, m_c2 );
441 double pStar0 = getPStar( m0, m_c1, m_c2 );
442 double width = width0 * pStar / pStar0 * m0 / m;
443 return width;
444}
445
446double EvtDTopipienu::getWidth1( double m, double m0, double m_c1, double m_c2, double width0,
447 double rBW ) {
448 double pStar = getPStar( m, m_c1, m_c2 );
449 double pStar0 = getPStar( m0, m_c1, m_c2 );
450 double F = getF1( m, m0, m_c1, m_c2, rBW );
451 double width = width0 * pStar / pStar0 * m0 / m * F * F;
452 return width;
453}
454
455double EvtDTopipienu::getWidth2( double m, double m0, double m_c1, double m_c2, double width0,
456 double rBW ) {
457 double pStar = getPStar( m, m_c1, m_c2 );
458 double pStar0 = getPStar( m0, m_c1, m_c2 );
459 double F = getF2( m, m0, m_c1, m_c2, rBW );
460 double width = width0 * pStar / pStar0 * m0 / m * F * F;
461 return width;
462}
463
464EvtComplex EvtDTopipienu::getCoef( double rho, double phi ) {
465 EvtComplex coef( rho * cos( phi ), rho * sin( phi ) );
466 return coef;
467}
468
469inline double EvtDTopipienu::getGx( double m0, double p0, double m_c1, double m_c2 ) {
470 double Gg = 0;
471 double MPI = 0.5 * ( m_c1 + m_c2 );
472 Gg = 3 * MPI * MPI * log( ( m0 + 2 * p0 ) / ( 2 * MPI ) ) / ( Pi * p0 * p0 ) +
473 m0 / ( 2 * Pi * p0 ) - MPI * MPI * m0 / ( Pi * p0 * p0 * p0 );
474 return Gg;
475}
476
477inline double EvtDTopipienu::getFx( double mr2, double sx, double p, double p0, double m_c1,
478 double m_c2 ) {
479 double Fx = 0;
480 Fx = mr2 / ( pow( p0, 3 ) ) *
481 ( p * p * ( getHx( sx, p, m_c1, m_c2 ) - getHx( mr2, p0, m_c1, m_c2 ) ) +
482 ( mr2 - sx ) * p0 * p0 * getdh( mr2, p0, m_c1, m_c2 ) );
483 return Fx;
484}
485
486inline double EvtDTopipienu::getHx( double sx, double p, double m_c1, double m_c2 ) {
487 double m = sqrt( sx );
488 double Hx = 0;
489 Hx = 2 * p * log( ( m + 2 * p ) / ( m_c1 + m_c2 ) ) / ( Pi * m );
490 return Hx;
491}
492
493inline double EvtDTopipienu::getdh( double mr2, double p0, double m_c1, double m_c2 ) {
494 double mass = sqrt( mr2 );
495 double dh =
496 getHx( mass, p0, m_c1, m_c2 ) * ( 1.0 / ( 8 * p0 * p0 ) - 1.0 / ( 2 * mass * mass ) ) +
497 1.0 / ( 2 * Pi * mass * mass );
498 return dh;
499}
500
501inline double EvtDTopipienu::getG1( double m2, double Mr ) {
502 double b1 = 1.302;
503 double b2 = 0.340;
504 double A = 2.426;
505 double Mr2 = Mr * Mr; // 0.953*0.953;
506 double gg1 = Mr * ( b1 + b2 * m2 ) * exp( -( m2 - Mr2 ) / A );
507 return gg1;
508}
509
510inline double EvtDTopipienu::getZ( double m2, double Mr2 ) {
511 double zz =
512 ( getRho( m2, mPi ) * log( ( 1.0 - getRho( m2, mPi ) ) / ( 1.0 + getRho( m2, mPi ) ) ) -
513 getRho( Mr2, mPi ) *
514 log( ( 1.0 - getRho( Mr2, mPi ) ) / ( 1.0 + getRho( Mr2, mPi ) ) ) ) /
515 Pi;
516
517 return zz;
518}
519
520inline double EvtDTopipienu::getRho( double m2, double mX ) {
521 double rho = 0.0;
522 if ( ( 1.0 - 4.0 * mX * mX / m2 ) > 0 ) rho = sqrt( 1.0 - 4.0 * mX * mX / m2 );
523 else rho = 0.0;
524 return rho;
525}
526
527inline EvtComplex EvtDTopipienu::getrho( double m2, double mX ) {
528 EvtComplex rho( 0.0, 0.0 );
529 EvtComplex ciR( 1.0, 0.0 );
530 EvtComplex ciM( 0.0, 1.0 );
531 if ( ( 1.0 - 4.0 * mX * mX / m2 ) > 0 ) rho = ciR * sqrt( 1.0 - 4.0 * mX * mX / m2 );
532 else rho = ciM * sqrt( 4.0 * mX * mX / m2 - 1.0 );
533 return rho;
534}
535inline double EvtDTopipienu::getWidthrho( double m, double m0, double width0, double p,
536 double p0 ) {
537 double widthRho = 0.0;
538 widthRho = width0 * pow( p / p0, 3 ) * ( m0 / m );
539 return widthRho;
540}
mass
double mass
Definition CosmicGenerator.cxx:128
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
pi2
character *LEPTONflag integer iresonances real pi2
Definition Ekhara/src/ekhara/common.h:3
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition Evt3Rank3C.cc:166
imag
double imag(const EvtComplex &c)
Definition EvtComplex.hh:207
abs2
double abs2(const EvtComplex &c)
Definition EvtComplex.hh:193
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition EvtDiracSpinor.cc:67
alpha
double alpha
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
EvtDTopipienu::getName
void getName(std::string &name)
Definition EvtDTopipienu.cc:33
EvtDTopipienu::decay
void decay(EvtParticle *p)
Definition EvtDTopipienu.cc:97
EvtDTopipienu::clone
EvtDecayBase * clone()
Definition EvtDTopipienu.cc:35
EvtDTopipienu::~EvtDTopipienu
virtual ~EvtDTopipienu()
Definition EvtDTopipienu.cc:31
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
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