EvtGen 2.2.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
Loading...
Searching...
No Matches
EvtDToKpienu.cpp
Go to the documentation of this file.
1
2/***********************************************************************
3* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4* *
5* This file is part of EvtGen. *
6* *
7* EvtGen is free software: you can redistribute it and/or modify *
8* it under the terms of the GNU General Public License as published by *
9* the Free Software Foundation, either version 3 of the License, or *
10* (at your option) any later version. *
11* *
12* EvtGen is distributed in the hope that it will be useful, *
13* but WITHOUT ANY WARRANTY; without even the implied warranty of *
14* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15* GNU General Public License for more details. *
16* *
17* You should have received a copy of the GNU General Public License *
18* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19***********************************************************************/
20
22
26#include "EvtGenBase/EvtPDL.hh"
29
30std::string EvtDToKpienu::getName() const
31{
32 return "DToKpienu";
33}
34
36{
37 return new EvtDToKpienu;
38}
39
41{
42 checkNArg( 0 );
43 checkNDaug( 4 );
47
48 EvtGenReport( EVTGEN_INFO, "EvtGen" )
49 << "EvtDToKpienu ==> Initialization !" << std::endl;
50 m_nAmps = 2;
51
52 m_rS = -11.57; // S-wave
53 m_rS1 = 0.08;
54 m_a_delta = 1.94;
55 m_b_delta = -0.81;
56 m_m0_1430_S = 1.425;
57 m_width0_1430_S = 0.270;
58 m_type[0] = 0;
59
60 m_mV = 1.81;
61 m_mA = 2.61;
62 m_V_0 = 1.411;
63 m_A1_0 = 1;
64 m_A2_0 = 0.788;
65
66 m_m0 = 0.8946; // P-wave K*
67 m_width0 = 0.04642;
68 m_rBW = 3.07;
69 m_rho = 1;
70 m_phi = 0;
71 m_type[1] = 1;
72
73 m_m0_1410 = 1.414; // P-wave K*(1410)
74 m_width0_1410 = 0.232;
75 m_rho_1410 = 0.1;
76 m_phi_1410 = 0.;
77 m_type[2] = 2;
78
79 m_TV_0 = 1; // D-wave K*2(1430)
80 m_T1_0 = 1;
81 m_T2_0 = 1;
82 m_m0_1430 = 1.4324;
83 m_width0_1430 = 0.109;
84 m_rho_1430 = 15;
85 m_phi_1430 = 0;
86 m_type[3] = 3;
87
88 m_mD = 1.86962;
89 m_mPi = 0.13957;
90 m_mK = 0.49368;
91 m_Pi = atan2( 0.0, -1.0 );
92 m_root2 = sqrt( 2. );
93 m_root2d3 = sqrt( 2. / 3 );
94 m_root1d2 = sqrt( 0.5 );
95 m_root3d2 = sqrt( 1.5 );
96}
97
99{
100 /*
101 * This piece of code could in principle be used to calculate maximum
102 * probablity on fly. But as it uses high number of points and model
103 * deals with single final state, we keep hardcoded number for now rather
104 * than adapting code to work here.
105
106 double maxprob = 0.0;
107 for(int ir=0;ir<=60000000;ir++){
108 p->initializePhaseSpace(getNDaug(),getDaugs());
109 EvtVector4R _K = p->getDaug(0)->getP4();
110 EvtVector4R _pi = p->getDaug(1)->getP4();
111 EvtVector4R _e = p->getDaug(2)->getP4();
112 EvtVector4R _nu = p->getDaug(3)->getP4();
113 int pid = EvtPDL::getStdHep(p->getDaug(0)->getId());
114 int charm;
115 if(pid == -321) charm = 1;
116 else charm = -1;
117 double m2, q2, cosV, cosL, chi;
118 KinVGen(_K, _pi, _e, _nu, charm, m2, q2, cosV, cosL, chi);
119 double _prob = calPDF(m2, q2, cosV, cosL, chi);
120 if(_prob>maxprob) {
121 maxprob=_prob;
122 EvtGenReport(EVTGEN_INFO,"EvtGen") << "Max PDF = " << ir << " charm= " << charm << " prob= "
123 << _prob << std::endl;
124 }
125 }
126 EvtGenReport(EVTGEN_INFO,"EvtGen") << "Max!!!!!!!!!!! " << maxprob<< std::endl;
127 */
128 setProbMax( 22750.0 );
129}
130
132{
134 EvtVector4R K = p->getDaug( 0 )->getP4();
135 EvtVector4R pi = p->getDaug( 1 )->getP4();
136 EvtVector4R e = p->getDaug( 2 )->getP4();
137 EvtVector4R nu = p->getDaug( 3 )->getP4();
138
139 int pid = EvtPDL::getStdHep( p->getDaug( 0 )->getId() );
140 int charm;
141 if ( pid == -321 ) {
142 charm = 1;
143 } else {
144 charm = -1;
145 }
146 double m2, q2, cosV, cosL, chi;
147 KinVGen( K, pi, e, nu, charm, m2, q2, cosV, cosL, chi );
148 double prob = calPDF( m2, q2, cosV, cosL, chi );
149 setProb( prob );
150 return;
151}
152
153void EvtDToKpienu::KinVGen( const EvtVector4R& vp4_K, const EvtVector4R& vp4_Pi,
154 const EvtVector4R& vp4_Lep, const EvtVector4R& vp4_Nu,
155 const int charm, double& m2, double& q2,
156 double& cosV, double& cosL, double& chi ) const
157{
158 EvtVector4R vp4_KPi = vp4_K + vp4_Pi;
159 EvtVector4R vp4_W = vp4_Lep + vp4_Nu;
160
161 m2 = vp4_KPi.mass2();
162 q2 = vp4_W.mass2();
163
164 EvtVector4R boost;
165 boost.set( vp4_KPi.get( 0 ), -vp4_KPi.get( 1 ), -vp4_KPi.get( 2 ),
166 -vp4_KPi.get( 3 ) );
167 EvtVector4R vp4_Kp = boostTo( vp4_K, boost );
168 cosV = vp4_Kp.dot( vp4_KPi ) / ( vp4_Kp.d3mag() * vp4_KPi.d3mag() );
169
170 boost.set( vp4_W.get( 0 ), -vp4_W.get( 1 ), -vp4_W.get( 2 ), -vp4_W.get( 3 ) );
171 EvtVector4R vp4_Lepp = boostTo( vp4_Lep, boost );
172 cosL = vp4_Lepp.dot( vp4_W ) / ( vp4_Lepp.d3mag() * vp4_W.d3mag() );
173
174 EvtVector4R V = vp4_KPi / vp4_KPi.d3mag();
175 EvtVector4R C = vp4_Kp.cross( V );
176 C /= C.d3mag();
177 EvtVector4R D = vp4_Lepp.cross( V );
178 D /= D.d3mag();
179 double sinx = C.cross( V ).dot( D );
180 double cosx = C.dot( D );
181 chi = sinx > 0 ? acos( cosx ) : -acos( cosx );
182 if ( charm == -1 )
183 chi = -chi;
184}
185
186double EvtDToKpienu::calPDF( const double m2, const double q2, const double cosV,
187 const double cosL, const double chi ) const
188{
189 double m = sqrt( m2 );
190 double q = sqrt( q2 );
191
192 // begin to calculate form factor
193 EvtComplex F10( 0.0, 0.0 );
194 EvtComplex F11( 0.0, 0.0 );
195 EvtComplex F21( 0.0, 0.0 );
196 EvtComplex F31( 0.0, 0.0 );
197 EvtComplex F12( 0.0, 0.0 );
198 EvtComplex F22( 0.0, 0.0 );
199 EvtComplex F32( 0.0, 0.0 );
200
201 EvtComplex f10( 0.0, 0.0 );
202 EvtComplex f11( 0.0, 0.0 );
203 EvtComplex f21( 0.0, 0.0 );
204 EvtComplex f31( 0.0, 0.0 );
205 EvtComplex f12( 0.0, 0.0 );
206 EvtComplex f22( 0.0, 0.0 );
207 EvtComplex f32( 0.0, 0.0 );
208 EvtComplex coef( 0.0, 0.0 );
209 double amplitude_temp, delta_temp;
210
211 for ( int index = 0; index < m_nAmps; index++ ) {
212 switch ( m_type[index] ) {
213 case 0: // calculate form factor of S wave
214 {
216 m_width0_1430_S, amplitude_temp, delta_temp, f10 );
217 F10 = F10 + f10;
218 break;
219 }
220 case 1: // calculate form factor of P wave (K*)
221 {
223 m_width0, m_rBW, amplitude_temp, delta_temp, f11,
224 f21, f31 );
225 coef = getCoef( m_rho, m_phi );
226 F11 = F11 + coef * f11;
227 F21 = F21 + coef * f21;
228 F31 = F31 + coef * f31;
229 break;
230 }
231 case 2: // calculate form factor of P wave (K*(1410))
232 {
234 m_width0_1410, m_rBW, amplitude_temp, delta_temp,
235 f11, f21, f31 );
236 coef = getCoef( m_rho_1410, m_phi_1410 );
237 F11 = F11 + coef * f11;
238 F21 = F21 + coef * f21;
239 F31 = F31 + coef * f31;
240 break;
241 }
242 case 3: // calculate form factor of D wave
243 {
245 m_width0_1430, m_rBW, amplitude_temp, delta_temp,
246 f12, f22, f32 );
247 coef = getCoef( m_rho_1430, m_phi_1430 );
248 F12 = F12 + coef * f12;
249 F22 = F22 + coef * f22;
250 F32 = F32 + coef * f32;
251 break;
252 }
253 default: {
254 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
255 << "No such form factor type!!!" << std::endl;
256 break;
257 }
258 }
259 }
260
261 // begin to calculate pdf value
262 double I, I1, I2, I3, I4, I5, I6, I7, I8, I9;
263
264 double cosV2 = cosV * cosV;
265 double sinV = sqrt( 1.0 - cosV2 );
266 double sinV2 = sinV * sinV;
267
268 EvtComplex F1 = F10 + F11 * cosV + F12 * ( 1.5 * cosV2 - 0.5 );
269 EvtComplex F2 = F21 * m_root1d2 + F22 * cosV * m_root3d2;
270 EvtComplex F3 = F31 * m_root1d2 + F32 * cosV * m_root3d2;
271
272 I1 = 0.25 * ( abs2( F1 ) + 1.5 * sinV2 * ( abs2( F2 ) + abs2( F3 ) ) );
273 I2 = -0.25 * ( abs2( F1 ) - 0.5 * sinV2 * ( abs2( F2 ) + abs2( F3 ) ) );
274 I3 = -0.25 * ( abs2( F2 ) - abs2( F3 ) ) * sinV2;
275 I4 = real( conj( F1 ) * F2 ) * sinV * 0.5;
276 I5 = real( conj( F1 ) * F3 ) * sinV;
277 I6 = real( conj( F2 ) * F3 ) * sinV2;
278 I7 = imag( conj( F2 ) * F1 ) * sinV;
279 I8 = imag( conj( F3 ) * F1 ) * sinV * 0.5;
280 I9 = imag( conj( F3 ) * F2 ) * sinV2 * ( -0.5 );
281
282 double coschi = cos( chi );
283 double sinchi = sin( chi );
284 double sin2chi = 2.0 * sinchi * coschi;
285 double cos2chi = 1.0 - 2.0 * sinchi * sinchi;
286
287 double sinL = sqrt( 1. - cosL * cosL );
288 double sinL2 = sinL * sinL;
289 double sin2L = 2.0 * sinL * cosL;
290 double cos2L = 1.0 - 2.0 * sinL2;
291
292 I = I1 + I2 * cos2L + I3 * sinL2 * cos2chi + I4 * sin2L * coschi +
293 I5 * sinL * coschi + I6 * cosL + I7 * sinL * sinchi +
294 I8 * sin2L * sinchi + I9 * sinL2 * sin2chi;
295 return I;
296}
297
298void EvtDToKpienu::ResonanceP( const double m, const double q, const double mV,
299 const double mA, const double V_0,
300 const double A1_0, const double A2_0,
301 const double m0, const double width0,
302 const double rBW, double& amplitude,
303 double& delta, EvtComplex& F11, EvtComplex& F21,
304 EvtComplex& F31 ) const
305{
306 double pKPi = getPStar( m_mD, m, q );
307 double mD2 = m_mD * m_mD;
308 double m2 = m * m;
309 double m02 = m0 * m0;
310 double q2 = q * q;
311 double mV2 = mV * mV;
312 double mA2 = mA * mA;
313 double summDm = m_mD + m;
314 double V = V_0 / ( 1.0 - q2 / ( mV2 ) );
315 double A1 = A1_0 / ( 1.0 - q2 / ( mA2 ) );
316 double A2 = A2_0 / ( 1.0 - q2 / ( mA2 ) );
317 double A = summDm * A1;
318 double B = 2.0 * m_mD * pKPi / summDm * V;
319
320 // construct the helicity form factor
321 double H0 = 0.5 / ( m * q ) *
322 ( ( mD2 - m2 - q2 ) * summDm * A1 -
323 4.0 * ( mD2 * pKPi * pKPi ) / summDm * A2 );
324 double Hp = A - B;
325 double Hm = A + B;
326
327 // calculate alpha
328 double B_Kstar = 2. / 3.; // B_Kstar = Br(Kstar(892)->k pi)
329 double pStar0 = getPStar( m0, m_mPi, m_mK );
330 double alpha = sqrt( 3. * m_Pi * B_Kstar / ( pStar0 * width0 ) );
331
332 // construct amplitudes of (non)resonance
333 double F = getF1( m, m0, m_mPi, m_mK, rBW );
334 double width = getWidth1( m, m0, m_mPi, m_mK, width0, rBW );
335
336 EvtComplex C( m0 * width0 * F, 0.0 );
337 double AA = m02 - m2;
338 double BB = -m0 * width;
339 EvtComplex amp = C / EvtComplex( AA, BB );
340 amplitude = abs( amp );
341 delta = atan2( imag( amp ), real( amp ) );
342
343 double alpham2 = alpha * 2.0;
344 F11 = amp * alpham2 * q * H0 * m_root2;
345 F21 = amp * alpham2 * q * ( Hp + Hm );
346 F31 = amp * alpham2 * q * ( Hp - Hm );
347}
348
349void EvtDToKpienu::NRS( const double m, const double q, const double rS,
350 const double rS1, const double a_delta,
351 const double b_delta, const double mA, const double m0,
352 const double width0, double& amplitude, double& delta,
353 EvtComplex& F10 ) const
354{
355 static const double tmp = ( m_mK + m_mPi ) * ( m_mK + m_mPi );
356
357 double m2 = m * m;
358 double q2 = q * q;
359 double mA2 = mA * mA;
360 double pKPi = getPStar( m_mD, m, q );
361 double m_K0_1430 = m0;
362 double width_K0_1430 = width0;
363 double m2_K0_1430 = m_K0_1430 * m_K0_1430;
364 double width = getWidth0( m, m_K0_1430, m_mPi, m_mK, width_K0_1430 );
365
366 // calculate modul of the amplitude
367 double x, Pm;
368 if ( m < m_K0_1430 ) {
369 x = sqrt( m2 / tmp - 1.0 );
370 Pm = 1.0 + rS1 * x;
371 } else {
372 x = sqrt( m2_K0_1430 / tmp - 1.0 );
373 Pm = 1.0 + rS1 * x;
374 Pm *= m_K0_1430 * width_K0_1430 /
375 sqrt( ( m2_K0_1430 - m2 ) * ( m2_K0_1430 - m2 ) +
376 m2_K0_1430 * width * width );
377 }
378
379 // calculate phase of the amplitude
380 double pStar = getPStar( m, m_mPi, m_mK );
381 double delta_bg = atan( 2. * a_delta * pStar /
382 ( 2. + a_delta * b_delta * pStar * pStar ) );
383 delta_bg = ( delta_bg > 0 ) ? delta_bg : ( delta_bg + m_Pi );
384
385 double delta_K0_1430 = atan( m_K0_1430 * width / ( m2_K0_1430 - m2 ) );
386 delta_K0_1430 = ( delta_K0_1430 > 0 ) ? delta_K0_1430
387 : ( delta_K0_1430 + m_Pi );
388 delta = delta_bg + delta_K0_1430;
389
390 EvtComplex ci( cos( delta ), sin( delta ) );
391 EvtComplex amp = ci * rS * Pm;
392 amplitude = rS * Pm;
393
394 F10 = amp * pKPi * m_mD / ( 1. - q2 / mA2 );
395}
396
397void EvtDToKpienu::ResonanceD( const double m, const double q, const double mV,
398 const double mA, const double TV_0,
399 const double T1_0, const double T2_0,
400 const double m0, const double width0,
401 const double rBW, double& amplitude,
402 double& delta, EvtComplex& F12, EvtComplex& F22,
403 EvtComplex& F32 ) const
404{
405 double pKPi = getPStar( m_mD, m, q );
406 double mD2 = m_mD * m_mD;
407 double m2 = m * m;
408 double m02 = m0 * m0;
409 double q2 = q * q;
410 double mV2 = mV * mV;
411 double mA2 = mA * mA;
412 double summDm = m_mD + m;
413 double TV = TV_0 / ( 1.0 - q2 / ( mV2 ) );
414 double T1 = T1_0 / ( 1.0 - q2 / ( mA2 ) );
415 double T2 = T2_0 / ( 1.0 - q2 / ( mA2 ) );
416
417 // construct amplitudes of (non)resonance
418 double F = getF2( m, m0, m_mPi, m_mK, rBW );
419 double width = getWidth2( m, m0, m_mPi, m_mK, width0, rBW );
420 EvtComplex C( m0 * width0 * F, 0.0 );
421 double AA = m02 - m2;
422 double BB = -m0 * width;
423
424 EvtComplex amp = C / EvtComplex( AA, BB );
425 amplitude = abs( amp );
426 delta = atan2( imag( amp ), real( amp ) );
427
428 F12 = amp * m_mD * pKPi / 3. *
429 ( ( mD2 - m2 - q2 ) * summDm * T1 - mD2 * pKPi * pKPi / summDm * T2 );
430 F22 = amp * m_root2d3 * m_mD * m * q * pKPi * summDm * T1;
431 F32 = amp * m_root2d3 * 2. * mD2 * m * q * pKPi * pKPi / summDm * TV;
432}
433
434double EvtDToKpienu::getPStar( const double m, const double m1,
435 const double m2 ) const
436{
437 double s = m * m;
438 double s1 = m1 * m1;
439 double s2 = m2 * m2;
440 double x = s + s1 - s2;
441 double t = 0.25 * x * x / s - s1;
442 double p;
443 if ( t > 0.0 ) {
444 p = sqrt( t );
445 } else {
446 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
447 << " Hello, pstar is less than 0.0" << std::endl;
448 p = 0.04;
449 }
450 return p;
451}
452
453double EvtDToKpienu::getF1( const double m, const double m0, const double m_c1,
454 const double m_c2, const double rBW ) const
455{
456 double pStar = getPStar( m, m_c1, m_c2 );
457 double pStar0 = getPStar( m0, m_c1, m_c2 );
458 double rBW2 = rBW * rBW;
459 double pStar2 = pStar * pStar;
460 double pStar02 = pStar0 * pStar0;
461 double B = 1. / sqrt( 1. + rBW2 * pStar2 );
462 double B0 = 1. / sqrt( 1. + rBW2 * pStar02 );
463 double F = pStar / pStar0 * B / B0;
464 return F;
465}
466
467double EvtDToKpienu::getF2( const double m, const double m0, const double m_c1,
468 const double m_c2, const double rBW ) const
469{
470 double pStar = getPStar( m, m_c1, m_c2 );
471 double pStar0 = getPStar( m0, m_c1, m_c2 );
472 double rBW2 = rBW * rBW;
473 double pStar2 = pStar * pStar;
474 double pStar02 = pStar0 * pStar0;
475 double B = 1. / sqrt( ( rBW2 * pStar2 - 3. ) * ( rBW2 * pStar2 - 3. ) +
476 9. * rBW2 * pStar2 );
477 double B0 = 1. / sqrt( ( rBW2 * pStar02 - 3. ) * ( rBW2 * pStar02 - 3. ) +
478 9. * rBW2 * pStar02 );
479 double F = pStar2 / pStar02 * B / B0;
480 return F;
481}
482
483double EvtDToKpienu::getWidth0( const double m, const double m0,
484 const double m_c1, const double m_c2,
485 const double width0 ) const
486{
487 double pStar = getPStar( m, m_c1, m_c2 );
488 double pStar0 = getPStar( m0, m_c1, m_c2 );
489 double width = width0 * pStar / pStar0 * m0 / m;
490 return width;
491}
492
493double EvtDToKpienu::getWidth1( const double m, const double m0,
494 const double m_c1, const double m_c2,
495 const double width0, const double rBW ) const
496{
497 double pStar = getPStar( m, m_c1, m_c2 );
498 double pStar0 = getPStar( m0, m_c1, m_c2 );
499 double F = getF1( m, m0, m_c1, m_c2, rBW );
500 double width = width0 * pStar / pStar0 * m0 / m * F * F;
501 return width;
502}
503
504double EvtDToKpienu::getWidth2( const double m, const double m0,
505 const double m_c1, const double m_c2,
506 const double width0, const double rBW ) const
507{
508 double pStar = getPStar( m, m_c1, m_c2 );
509 double pStar0 = getPStar( m0, m_c1, m_c2 );
510 double F = getF2( m, m0, m_c1, m_c2, rBW );
511 double width = width0 * pStar / pStar0 * m0 / m * F * F;
512 return width;
513}
514
515EvtComplex EvtDToKpienu::getCoef( const double rho, const double phi ) const
516{
517 EvtComplex coef( rho * cos( phi ), rho * sin( phi ) );
518 return coef;
519}
const float pi
const EvtComplex I
EvtComplex conj(const EvtComplex &c)
double imag(const EvtComplex &c)
double real(const EvtComplex &c)
double abs2(const EvtComplex &c)
double abs(const EvtComplex &c)
EvtRaritaSchwinger boostTo(const EvtRaritaSchwinger &rs, const EvtVector4R p4)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=nullptr)
Definition EvtReport.cpp:32
@ EVTGEN_INFO
Definition EvtReport.hh:52
@ EVTGEN_ERROR
Definition EvtReport.hh:49
double m_width0_1430_S
double getWidth0(const double m, const double m0, const double m_c1, const double m_c2, const double width0) const
double m_rho_1410
double getPStar(const double m, const double m1, const double m2) const
void initProbMax() override
void KinVGen(const EvtVector4R &vp4_K, const EvtVector4R &vp4_Pi, const EvtVector4R &vp4_Lep, const EvtVector4R &vp4_Nu, const int charm, double &m2, double &q2, double &cosV, double &cosL, double &chi) const
void ResonanceD(const double m, const double q, const double mV, const double mA, const double TV_0, const double T1_0, const double T2_0, const double m0, const double width0, const double rBW, double &amplitude, double &delta, EvtComplex &F12, EvtComplex &F22, EvtComplex &F32) const
double m_width0_1410
double getWidth2(const double m, const double m0, const double m_c1, const double m_c2, const double width0, const double rBW) const
void init() override
EvtDecayBase * clone() const override
double getF2(const double m, const double m0, const double m_c1, const double m_c2, const double rBW) const
double m_m0_1430_S
double getWidth1(const double m, const double m0, const double m_c1, const double m_c2, const double width0, const double rBW) const
void ResonanceP(const double m, const double q, const double mV, const double mA, const double V_0, const double A1_0, const double A2_0, const double m0, const double width0, const double rBW, double &amplitude, double &delta, EvtComplex &F11, EvtComplex &F21, EvtComplex &F31) const
void decay(EvtParticle *p) override
void NRS(const double m, const double q, const double rS, const double rS1, const double a_delta, const double b_delta, const double mA, const double m0, const double width0, double &amplitude, double &delta, EvtComplex &F10) const
std::array< int, 5 > m_type
double getF1(const double m, const double m0, const double m_c1, const double m_c2, const double rBW) const
double m_width0_1430
double calPDF(const double m2, const double q2, const double cosV, const double cosL, const double chi) const
EvtComplex getCoef(const double rho, const double phi) const
std::string getName() const override
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
EvtDecayBase()=default
int getNDaug() const
void checkSpinParent(EvtSpinType::spintype sp)
void setProbMax(double prbmx)
void checkNDaug(int d1, int d2=-1)
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
const EvtId * getDaugs() const
void setProb(double prob)
static int getStdHep(EvtId id)
Definition EvtPDL.cpp:356
double initializePhaseSpace(size_t numdaughter, const EvtId *daughters, bool forceResetMasses=false, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
EvtId getId() const
const EvtVector4R & getP4() const
EvtParticle * getDaug(const int i)
double dot(const EvtVector4R &v2) const
double get(int i) const
double d3mag() const
EvtVector4R cross(const EvtVector4R &v2) const
double mass2() const
void set(int i, double d)