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