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
EvtDalitzReso.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"
30
31#include <assert.h>
32#include <cmath>
33#include <iostream>
34#include <stdlib.h>
35
36#define PRECISION ( 1.e-3 )
37
40
41// single Breit-Wigner
43 EvtSpinType::spintype spin, double m0, double g0,
44 NumType typeN, double f_b, double f_d ) :
45 m_dp( dp ),
46 m_pairAng( pairAng ),
47 m_pairRes( pairRes ),
48 m_spin( spin ),
49 m_typeN( typeN ),
50 m_m0( m0 ),
51 m_g0( g0 ),
52 m_massFirst( dp.m( first( pairRes ) ) ),
53 m_massSecond( dp.m( second( pairRes ) ) ),
54 m_m0_mix( -1. ),
55 m_g0_mix( 0. ),
56 m_delta_mix( 0. ),
57 m_amp_mix( 0., 0. ),
58 m_g1( -1. ),
59 m_g2( -1. ),
61 m_f_b( f_b ),
62 m_f_d( f_d ),
63 m_kmatrix_index( -1 ),
64 m_fr12prod( 0., 0. ),
65 m_fr13prod( 0., 0. ),
66 m_fr14prod( 0., 0. ),
67 m_fr15prod( 0., 0. ),
68 m_s0prod( 0. ),
69 m_a( 0. ),
70 m_r( 0. ),
71 m_Blass( 0. ),
72 m_phiB( 0. ),
73 m_R( 0. ),
74 m_phiR( 0. ),
75 m_cutoff( -1. ),
76 m_scaleByMOverQ( false ),
77 m_alpha( 0. )
78{
79 assert( m_typeN != K_MATRIX && m_typeN != K_MATRIX_I &&
80 m_typeN != K_MATRIX_II ); // single BW cannot be K-matrix
81}
82
83// Breit-Wigner with electromagnetic mass mixing
85 EvtSpinType::spintype spin, double m0, double g0,
86 NumType typeN, double m0_mix, double g0_mix,
87 double delta_mix, EvtComplex amp_mix ) :
88 m_dp( dp ),
89 m_pairAng( pairAng ),
90 m_pairRes( pairRes ),
91 m_spin( spin ),
92 m_typeN( typeN ),
93 m_m0( m0 ),
94 m_g0( g0 ),
95 m_massFirst( dp.m( first( pairRes ) ) ),
96 m_massSecond( dp.m( second( pairRes ) ) ),
97 m_m0_mix( m0_mix ),
98 m_g0_mix( g0_mix ),
99 m_delta_mix( delta_mix ),
100 m_amp_mix( amp_mix ),
101 m_g1( -1. ),
102 m_g2( -1. ),
104 m_f_b( 0.0 ),
105 m_f_d( 1.5 ),
106 m_kmatrix_index( -1 ),
107 m_fr12prod( 0., 0. ),
108 m_fr13prod( 0., 0. ),
109 m_fr14prod( 0., 0. ),
110 m_fr15prod( 0., 0. ),
111 m_s0prod( 0. ),
112 m_a( 0. ),
113 m_r( 0. ),
114 m_Blass( 0. ),
115 m_phiB( 0. ),
116 m_R( 0. ),
117 m_phiR( 0. ),
118 m_cutoff( -1. ),
119 m_scaleByMOverQ( false ),
120 m_alpha( 0. )
121{
122 // single BW (with electromagnetic mixing) cannot be K-matrix
123 assert( m_typeN != K_MATRIX && m_typeN != K_MATRIX_I &&
124 m_typeN != K_MATRIX_II );
125}
126
127// coupled Breit-Wigner
129 Pair pairRes, EvtSpinType::spintype spin,
130 double m0, NumType typeN, double g1, double g2,
131 CouplingType coupling2 ) :
132 m_dp( dp ),
133 m_pairAng( pairAng ),
134 m_pairRes( pairRes ),
135 m_spin( spin ),
136 m_typeN( typeN ),
137 m_m0( m0 ),
138 m_g0( -1. ),
139 m_massFirst( dp.m( first( pairRes ) ) ),
140 m_massSecond( dp.m( second( pairRes ) ) ),
141 m_m0_mix( -1. ),
142 m_g0_mix( 0. ),
143 m_delta_mix( 0. ),
144 m_amp_mix( 0., 0. ),
145 m_g1( g1 ),
146 m_g2( g2 ),
147 m_coupling2( coupling2 ),
148 m_f_b( 0.0 ),
149 m_f_d( 1.5 ),
150 m_kmatrix_index( -1 ),
151 m_fr12prod( 0., 0. ),
152 m_fr13prod( 0., 0. ),
153 m_fr14prod( 0., 0. ),
154 m_fr15prod( 0., 0. ),
155 m_s0prod( 0. ),
156 m_a( 0. ),
157 m_r( 0. ),
158 m_Blass( 0. ),
159 m_phiB( 0. ),
160 m_R( 0. ),
161 m_phiR( 0. ),
162 m_cutoff( -1. ),
163 m_scaleByMOverQ( false ),
164 m_alpha( 0. )
165{
166 assert( m_coupling2 != Undefined );
167 assert( m_typeN != K_MATRIX && m_typeN != K_MATRIX_I &&
168 m_typeN != K_MATRIX_II ); // coupled BW cannot be K-matrix
169 assert( m_typeN != LASS ); // coupled BW cannot be LASS
170 assert( m_typeN != NBW ); // for coupled BW, only relativistic BW
171}
172
173// K-Matrix (A&S)
175 std::string nameIndex, NumType typeN,
176 EvtComplex fr12prod, EvtComplex fr13prod,
177 EvtComplex fr14prod, EvtComplex fr15prod,
178 double s0prod ) :
179 m_dp( dp ),
180 m_pairRes( pairRes ),
181 m_typeN( typeN ),
182 m_m0( 0. ),
183 m_g0( 0. ),
184 m_massFirst( dp.m( first( pairRes ) ) ),
185 m_massSecond( dp.m( second( pairRes ) ) ),
186 m_m0_mix( -1. ),
187 m_g0_mix( 0. ),
188 m_delta_mix( 0. ),
189 m_amp_mix( 0., 0. ),
190 m_g1( -1. ),
191 m_g2( -1. ),
193 m_f_b( 0. ),
194 m_f_d( 0. ),
195 m_kmatrix_index( -1 ),
196 m_fr12prod( fr12prod ),
197 m_fr13prod( fr13prod ),
198 m_fr14prod( fr14prod ),
199 m_fr15prod( fr15prod ),
200 m_s0prod( s0prod ),
201 m_a( 0. ),
202 m_r( 0. ),
203 m_Blass( 0. ),
204 m_phiB( 0. ),
205 m_R( 0. ),
206 m_phiR( 0. ),
207 m_cutoff( -1. ),
208 m_scaleByMOverQ( false ),
209 m_alpha( 0. )
210{
211 assert( m_typeN == K_MATRIX || m_typeN == K_MATRIX_I ||
212 m_typeN == K_MATRIX_II );
214 if ( nameIndex == "Pole1" )
215 m_kmatrix_index = 1;
216 else if ( nameIndex == "Pole2" )
217 m_kmatrix_index = 2;
218 else if ( nameIndex == "Pole3" )
219 m_kmatrix_index = 3;
220 else if ( nameIndex == "Pole4" )
221 m_kmatrix_index = 4;
222 else if ( nameIndex == "Pole5" )
223 m_kmatrix_index = 5;
224 else if ( nameIndex == "f11prod" )
225 m_kmatrix_index = 6;
226 else
227 assert( 0 );
228}
229
230// LASS parameterization
231EvtDalitzReso::EvtDalitzReso( const EvtDalitzPlot& dp, Pair pairRes, double m0,
232 double g0, double a, double r, double B,
233 double phiB, double R, double phiR, double cutoff,
234 bool scaleByMOverQ ) :
235 m_dp( dp ),
236 m_pairRes( pairRes ),
237 m_typeN( LASS ),
238 m_m0( m0 ),
239 m_g0( g0 ),
240 m_massFirst( dp.m( first( pairRes ) ) ),
241 m_massSecond( dp.m( second( pairRes ) ) ),
242 m_m0_mix( -1. ),
243 m_g0_mix( 0. ),
244 m_delta_mix( 0. ),
245 m_amp_mix( 0., 0. ),
246 m_g1( -1. ),
247 m_g2( -1. ),
249 m_f_b( 0.0 ),
250 m_f_d( 1.5 ),
251 m_kmatrix_index( -1 ),
252 m_fr12prod( 0., 0. ),
253 m_fr13prod( 0., 0. ),
254 m_fr14prod( 0., 0. ),
255 m_fr15prod( 0., 0. ),
256 m_s0prod( 0. ),
257 m_a( a ),
258 m_r( r ),
259 m_Blass( B ),
260 m_phiB( phiB ),
261 m_R( R ),
262 m_phiR( phiR ),
263 m_cutoff( cutoff ),
264 m_scaleByMOverQ( scaleByMOverQ ),
265 m_alpha( 0. )
266{
268}
269
270//Flatte
272 double m0 ) :
273 m_dp( dp ),
274 m_pairRes( pairRes ),
275 m_typeN( FLATTE ),
276 m_m0( m0 ),
277 m_g0( 0. ),
278 m_massFirst( dp.m( first( pairRes ) ) ),
279 m_massSecond( dp.m( second( pairRes ) ) ),
280 m_m0_mix( -1. ),
281 m_g0_mix( 0. ),
282 m_delta_mix( 0. ),
283 m_amp_mix( 0., 0. ),
284 m_g1( -1. ),
285 m_g2( -1. ),
287 m_f_b( 0. ),
288 m_f_d( 0. ),
289 m_kmatrix_index( -1 ),
290 m_fr12prod( 0., 0. ),
291 m_fr13prod( 0., 0. ),
292 m_fr14prod( 0., 0. ),
293 m_fr15prod( 0., 0. ),
294 m_s0prod( 0. ),
295 m_a( 0. ),
296 m_r( 0. ),
297 m_Blass( 0. ),
298 m_phiB( 0. ),
299 m_R( 0. ),
300 m_phiR( 0. ),
301 m_cutoff( -1. ),
302 m_scaleByMOverQ( false ),
303 m_alpha( 0. )
304{
306}
307
309{
310 if ( m_typeN == NON_RES )
311 return EvtComplex( 1.0, 0.0 );
312
313 const double s = x.q( m_pairRes );
314 const double m = sqrt( s );
315
316 if ( m_typeN == NON_RES_LIN )
317 return s;
318
319 if ( m_typeN == NON_RES_EXP )
320 return exp( -m_alpha * s );
321
322 // do use always hash table (speed up fitting)
324 return Fvector( s, m_kmatrix_index );
325
326 if ( m_typeN == FLATTE )
327 return flatte( s );
328
330 vd.set_f( m_f_d );
331
333
334 if ( m_typeN == LASS )
335 return lass( kd, vd );
336
337 EvtComplex amp( 1.0, 0.0 );
338
340 x.bigM(), m_spin );
341 vb.set_f( m_f_b );
342
343 const EvtTwoBodyKine kb( m, x.m( EvtCyclic3::other( m_pairRes ) ), x.bigM() );
344
345 EvtComplex prop( 0, 0 );
346 if ( m_typeN == NBW ) {
347 prop = propBreitWigner( m_m0, m_g0, m );
348 } else if ( m_typeN == GAUSS_CLEO || m_typeN == GAUSS_CLEO_ZEMACH ) {
349 prop = propGauss( m_m0, m_g0, m );
350 } else {
351 if ( m_coupling2 == Undefined ) {
352 // single BW
353 const double g = ( m_g0 <= 0. || vd.pD() <= 0. )
354 ? -m_g0
355 : m_g0 * vd.widthFactor( kd ); // running width
356 if ( m_typeN == GS_CLEO || m_typeN == GS_CLEO_ZEMACH ) {
357 // Gounaris-Sakurai (GS)
358 prop = propGounarisSakurai( m_m0, fabs( m_g0 ), vd.pD(), m, g,
359 kd.p() );
360 } else {
361 // standard relativistic BW
362 prop = propBreitWignerRel( m_m0, g, m );
363 }
364 } else {
365 // coupled width BW
366 EvtComplex G1, G2;
367 switch ( m_coupling2 ) {
368 case PicPic: {
369 G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
370 static const double mPic = EvtPDL::getMass(
371 EvtPDL::getId( "pi+" ) );
372 G2 = m_g2 * m_g2 * psFactor( mPic, mPic, m );
373 break;
374 }
375 case PizPiz: {
376 G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
377 static const double mPiz = EvtPDL::getMass(
378 EvtPDL::getId( "pi0" ) );
379 G2 = m_g2 * m_g2 * psFactor( mPiz, mPiz, m );
380 break;
381 }
382 case PiPi: {
383 G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
384 static const double mPic = EvtPDL::getMass(
385 EvtPDL::getId( "pi+" ) );
386 static const double mPiz = EvtPDL::getMass(
387 EvtPDL::getId( "pi0" ) );
388 G2 = m_g2 * m_g2 * psFactor( mPic, mPic, mPiz, mPiz, m );
389 break;
390 }
391 case KcKc: {
392 G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
393 static const double mKc = EvtPDL::getMass(
394 EvtPDL::getId( "K+" ) );
395 G2 = m_g2 * m_g2 * psFactor( mKc, mKc, m );
396 break;
397 }
398 case KzKz: {
399 G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
400 static const double mKz = EvtPDL::getMass(
401 EvtPDL::getId( "K0" ) );
402 G2 = m_g2 * m_g2 * psFactor( mKz, mKz, m );
403 break;
404 }
405 case KK: {
406 G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
407 static const double mKc = EvtPDL::getMass(
408 EvtPDL::getId( "K+" ) );
409 static const double mKz = EvtPDL::getMass(
410 EvtPDL::getId( "K0" ) );
411 G2 = m_g2 * m_g2 * psFactor( mKc, mKc, mKz, mKz, m );
412 break;
413 }
414 case EtaPic: {
415 G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
416 static const double mEta = EvtPDL::getMass(
417 EvtPDL::getId( "eta" ) );
418 static const double mPic = EvtPDL::getMass(
419 EvtPDL::getId( "pi+" ) );
420 G2 = m_g2 * m_g2 * psFactor( mEta, mPic, m );
421 break;
422 }
423 case EtaPiz: {
424 G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
425 static const double mEta = EvtPDL::getMass(
426 EvtPDL::getId( "eta" ) );
427 static const double mPiz = EvtPDL::getMass(
428 EvtPDL::getId( "pi0" ) );
429 G2 = m_g2 * m_g2 * psFactor( mEta, mPiz, m );
430 break;
431 }
432 case PicPicKK: {
433 static const double mPic = EvtPDL::getMass(
434 EvtPDL::getId( "pi+" ) );
435 G1 = m_g1 * psFactor( mPic, mPic, m );
436 static const double mKc = EvtPDL::getMass(
437 EvtPDL::getId( "K+" ) );
438 static const double mKz = EvtPDL::getMass(
439 EvtPDL::getId( "K0" ) );
440 G2 = m_g2 * psFactor( mKc, mKc, mKz, mKz, m );
441 break;
442 }
443 default:
444 std::cout
445 << "EvtDalitzReso:evaluate(): PANIC, wrong coupling2 state."
446 << std::endl;
447 assert( 0 );
448 break;
449 }
450 // calculate standard couple BW propagator
451 if ( m_coupling2 != WA76 )
452 prop = m_g1 * propBreitWignerRelCoupled( m_m0, G1, G2, m );
453 }
454 }
455 amp *= prop;
456
457 // Compute form-factors (Blatt-Weisskopf penetration factor)
458 amp *= vb.formFactor( kb );
459 amp *= vd.formFactor( kd );
460
461 // Compute numerator (angular distribution)
462 amp *= numerator( x, vb, vd, kb, kd );
463
464 // Compute electromagnetic mass mixing factor
465 if ( m_m0_mix > 0. ) {
466 EvtComplex prop_mix;
467 if ( m_typeN == NBW ) {
468 prop_mix = propBreitWigner( m_m0_mix, m_g0_mix, m );
469 } else {
470 assert( m_g1 < 0. ); // running width only
471 const double g_mix = m_g0_mix * vd.widthFactor( kd );
472 prop_mix = propBreitWignerRel( m_m0_mix, g_mix, m );
473 }
474 amp *= mixFactor( prop, prop_mix );
475 }
476
477 return amp;
478}
479
480EvtComplex EvtDalitzReso::psFactor( const double ma, const double mb,
481 const double m ) const
482{
483 if ( m > ( ma + mb ) ) {
484 const EvtTwoBodyKine kd( ma, mb, m );
485 return EvtComplex( 0, 2 * kd.p() / m );
486 } else {
487 // analytical continuation
488 const double s = m * m;
489 const double phaseFactor_analyticalCont =
490 -0.5 * ( sqrt( 4 * ma * ma / s - 1 ) + sqrt( 4 * mb * mb / s - 1 ) );
491 return EvtComplex( phaseFactor_analyticalCont, 0 );
492 }
493}
494
495EvtComplex EvtDalitzReso::psFactor( const double ma1, const double mb1,
496 const double ma2, const double mb2,
497 const double m ) const
498{
499 return 0.5 * ( psFactor( ma1, mb1, m ) + psFactor( ma2, mb2, m ) );
500}
501
502EvtComplex EvtDalitzReso::propGauss( const double m0, const double s0,
503 const double m ) const
504{
505 // Gaussian
506 const double gauss = 1. / sqrt( EvtConst::twoPi ) / s0 *
507 exp( -( m - m0 ) * ( m - m0 ) / 2. / ( s0 * s0 ) );
508 return EvtComplex( gauss, 0. );
509}
510
511EvtComplex EvtDalitzReso::propBreitWigner( const double m0, const double g0,
512 const double m ) const
513{
514 // non-relativistic BW
515 return sqrt( g0 / EvtConst::twoPi ) / ( m - m0 - EvtComplex( 0.0, g0 / 2. ) );
516}
517
518EvtComplex EvtDalitzReso::propBreitWignerRel( const double m0, const double g0,
519 const double m ) const
520{
521 // relativistic BW with real width
522 return 1. / ( m0 * m0 - m * m - EvtComplex( 0., m0 * g0 ) );
523}
524
526 const EvtComplex& g0,
527 const double m ) const
528{
529 // relativistic BW with complex width
530 return 1. / ( m0 * m0 - m * m - EvtComplex( 0., m0 ) * g0 );
531}
532
534 const EvtComplex& g1,
535 const EvtComplex& g2,
536 const double m ) const
537{
538 // relativistic coupled BW
539 return 1. / ( m0 * m0 - m * m - ( g1 + g2 ) );
540}
541
542EvtComplex EvtDalitzReso::propGounarisSakurai( const double m0, const double g0,
543 const double k0, const double m,
544 const double g,
545 const double k ) const
546{
547 // Gounaris-Sakurai parameterization of pi+pi- P wave. PRD, Vol61, 112002. PRL, Vol21, 244.
548 // Expressions taken from BAD637v4, after fixing the imaginary part of the BW denominator: i M_R Gamma_R(s) --> i sqrt(s) Gamma_R(s)
549 return ( 1. + GS_d( m0, k0 ) * g0 / m0 ) /
550 ( m0 * m0 - m * m - EvtComplex( 0., m * g ) +
551 GS_f( m0, g0, k0, m, k ) );
552}
553
554inline double EvtDalitzReso::GS_f( const double m0, const double g0,
555 const double k0, const double m,
556 const double k ) const
557{
558 // m: sqrt(s)
559 // m0: nominal resonance mass
560 // k: momentum of pion in resonance rest frame (at m)
561 // k0: momentum of pion in resonance rest frame (at nominal resonance mass)
562 return g0 * m0 * m0 / ( k0 * k0 * k0 ) *
563 ( k * k * ( GS_h( m, k ) - GS_h( m0, k0 ) ) +
564 ( m0 * m0 - m * m ) * k0 * k0 * GS_dhods( m0, k0 ) );
565}
566
567inline double EvtDalitzReso::GS_h( const double m, const double k ) const
568{
569 return 2. / EvtConst::pi * k / m *
570 log( ( m + 2. * k ) / ( 2. * m_massFirst ) );
571}
572
573inline double EvtDalitzReso::GS_dhods( const double m0, const double k0 ) const
574{
575 return GS_h( m0, k0 ) * ( 0.125 / ( k0 * k0 ) - 0.5 / ( m0 * m0 ) ) +
576 0.5 / ( EvtConst::pi * m0 * m0 );
577}
578
579inline double EvtDalitzReso::GS_d( const double m0, const double k0 ) const
580{
581 return 3. / EvtConst::pi * m_massFirst * m_massFirst / ( k0 * k0 ) *
582 log( ( m0 + 2. * k0 ) / ( 2. * m_massFirst ) ) +
583 m0 / ( 2. * EvtConst::pi * k0 ) -
584 m_massFirst * m_massFirst * m0 / ( EvtConst::pi * k0 * k0 * k0 );
585}
586
588 const EvtTwoBodyVertex& vb,
589 const EvtTwoBodyVertex& vd,
590 const EvtTwoBodyKine& kb,
591 const EvtTwoBodyKine& kd ) const
592{
593 EvtComplex ret( 0., 0. );
594
595 // Non-relativistic Breit-Wigner
596 if ( NBW == m_typeN ) {
597 ret = angDep( x );
598 }
599
600 // Standard relativistic Zemach propagator
601 else if ( RBW_ZEMACH == m_typeN ) {
602 ret = vd.phaseSpaceFactor( kd, EvtTwoBodyKine::AB ) * angDep( x );
603 }
604
605 // Standard relativistic Zemach propagator
606 else if ( RBW_ZEMACH2 == m_typeN ) {
607 ret = vd.phaseSpaceFactor( kd, EvtTwoBodyKine::AB ) *
609 if ( m_spin == EvtSpinType::VECTOR ) {
610 ret *= -4.;
611 } else if ( m_spin == EvtSpinType::TENSOR ) {
612 ret *= 16. / 3.;
613 } else if ( m_spin != EvtSpinType::SCALAR )
614 assert( 0 );
615 }
616
617 // Kuehn-Santamaria normalization:
618 else if ( RBW_KUEHN == m_typeN ) {
619 ret = m_m0 * m_m0 * angDep( x );
620 }
621
622 // CLEO amplitude
623 else if ( ( RBW_CLEO == m_typeN ) || ( GS_CLEO == m_typeN ) ||
625 ( GAUSS_CLEO == m_typeN ) || ( GAUSS_CLEO_ZEMACH == m_typeN ) ) {
626 const Index iA = other( m_pairAng ); // A = other(BC)
627 const Index iB = common( m_pairRes, m_pairAng ); // B = common(AB,BC)
628 const Index iC = other( m_pairRes ); // C = other(AB)
629
630 const double M = x.bigM();
631 const double mA = x.m( iA );
632 const double mB = x.m( iB );
633 const double mC = x.m( iC );
634 const double qAB = x.q( combine( iA, iB ) );
635 const double qBC = x.q( combine( iB, iC ) );
636 const double qCA = x.q( combine( iC, iA ) );
637
638 const double M2 = M * M;
639 const double m02 = ( ( RBW_CLEO_ZEMACH == m_typeN ) ||
640 ( GS_CLEO_ZEMACH == m_typeN ) ||
642 ? qAB
643 : m_m0 * m_m0;
644 const double mA2 = mA * mA;
645 const double mB2 = mB * mB;
646 const double mC2 = mC * mC;
647
649 ret = EvtComplex( 1., 0. );
650 else if ( m_spin == EvtSpinType::VECTOR ) {
651 ret = qCA - qBC + ( M2 - mC2 ) * ( mB2 - mA2 ) / m02;
652 } else if ( m_spin == EvtSpinType::TENSOR ) {
653 const double x1 = qBC - qCA + ( M2 - mC2 ) * ( mA2 - mB2 ) / m02;
654 const double x2 = M2 - mC2;
655 const double x3 = qAB - 2 * M2 - 2 * mC2 + x2 * x2 / m02;
656 const double x4 = mA2 - mB2;
657 const double x5 = qAB - 2 * mB2 - 2 * mA2 + x4 * x4 / m02;
658 ret = x1 * x1 - x3 * x5 / 3.;
659 } else
660 assert( 0 );
661 }
662
663 return ret;
664}
665
666double EvtDalitzReso::angDep( const EvtDalitzPoint& x ) const
667{
668 // Angular dependece for factorizable amplitudes
669 // unphysical cosines indicate we are in big trouble
670 const double cosTh = x.cosTh(
671 m_pairAng, m_pairRes ); // angle between common(reso,ang) and other(reso)
672 if ( fabs( cosTh ) > 1. ) {
673 EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "cosTh " << cosTh << std::endl;
674 assert( 0 );
675 }
676
677 // in units of half-spin
678 return EvtdFunction::d( EvtSpinType::getSpin2( m_spin ), 2 * 0, 2 * 0,
679 acos( cosTh ) );
680}
681
683 const EvtComplex& prop_mix ) const
684{
685 const double Delta = m_delta_mix * ( m_m0 + m_m0_mix );
686 return 1 / ( 1 - Delta * Delta * prop * prop_mix ) *
687 ( 1 + m_amp_mix * Delta * prop_mix );
688}
689
690EvtComplex EvtDalitzReso::Fvector( const double s, const int index ) const
691{
692 assert( index >= 1 && index <= 6 );
693
694 //Define the complex coupling constant
695 //The convection is as follow
696 //i=0 --> pi+ pi-
697 //i=1 --> KK
698 //i=2 --> 4pi
699 //i=3 --> eta eta
700 //i=4 --> eta eta'
701 //The first index is the resonace-pole index
702
703 double g[5][5]; // Coupling constants. The first index is the pole index. The second index is the decay channel
704 double ma[5]; // Pole masses. The unit is in GeV
705
706 const int solution = ( m_typeN == K_MATRIX )
707 ? 3
708 : ( ( m_typeN == K_MATRIX_I )
709 ? 1
710 : ( ( m_typeN == K_MATRIX_II ) ? 2 : 0 ) );
711 if ( solution == 0 ) {
712 std::cout << "EvtDalitzReso::Fvector() error. Kmatrix solution incorrectly chosen ! "
713 << std::endl;
714 abort();
715 }
716
717 if ( solution == 3 ) {
718 // coupling constants
719 //pi+pi- channel
720 g[0][0] = 0.22889;
721 g[1][0] = 0.94128;
722 g[2][0] = 0.36856;
723 g[3][0] = 0.33650;
724 g[4][0] = 0.18171;
725 //K+K- channel
726 g[0][1] = -0.55377;
727 g[1][1] = 0.55095;
728 g[2][1] = 0.23888;
729 g[3][1] = 0.40907;
730 g[4][1] = -0.17558;
731 //4pi channel
732 g[0][2] = 0;
733 g[1][2] = 0;
734 g[2][2] = 0.55639;
735 g[3][2] = 0.85679;
736 g[4][2] = -0.79658;
737 //eta eta channel
738 g[0][3] = -0.39899;
739 g[1][3] = 0.39065;
740 g[2][3] = 0.18340;
741 g[3][3] = 0.19906;
742 g[4][3] = -0.00355;
743 //eta eta' channel
744 g[0][4] = -0.34639;
745 g[1][4] = 0.31503;
746 g[2][4] = 0.18681;
747 g[3][4] = -0.00984;
748 g[4][4] = 0.22358;
749
750 // Pole masses
751 ma[0] = 0.651;
752 ma[1] = 1.20360;
753 ma[2] = 1.55817;
754 ma[3] = 1.21000;
755 ma[4] = 1.82206;
756
757 } else if ( solution == 1 ) { // solnI.txt
758
759 // coupling constants
760 //pi+pi- channel
761 g[0][0] = 0.31896;
762 g[1][0] = 0.85963;
763 g[2][0] = 0.47993;
764 g[3][0] = 0.45121;
765 g[4][0] = 0.39391;
766 //K+K- channel
767 g[0][1] = -0.49998;
768 g[1][1] = 0.52402;
769 g[2][1] = 0.40254;
770 g[3][1] = 0.42769;
771 g[4][1] = -0.30860;
772 //4pi channel
773 g[0][2] = 0;
774 g[1][2] = 0;
775 g[2][2] = 1.0;
776 g[3][2] = 1.15088;
777 g[4][2] = 0.33999;
778 //eta eta channel
779 g[0][3] = -0.21554;
780 g[1][3] = 0.38093;
781 g[2][3] = 0.21811;
782 g[3][3] = 0.22925;
783 g[4][3] = 0.06919;
784 //eta eta' channel
785 g[0][4] = -0.18294;
786 g[1][4] = 0.23788;
787 g[2][4] = 0.05454;
788 g[3][4] = 0.06444;
789 g[4][4] = 0.32620;
790
791 // Pole masses
792 ma[0] = 0.7369;
793 ma[1] = 1.24347;
794 ma[2] = 1.62681;
795 ma[3] = 1.21900;
796 ma[4] = 1.74932;
797
798 } else if ( solution == 2 ) { // solnIIa.txt
799
800 // coupling constants
801 //pi+pi- channel
802 g[0][0] = 0.26014;
803 g[1][0] = 0.95289;
804 g[2][0] = 0.46244;
805 g[3][0] = 0.41848;
806 g[4][0] = 0.01804;
807 //K+K- channel
808 g[0][1] = -0.57849;
809 g[1][1] = 0.55887;
810 g[2][1] = 0.31712;
811 g[3][1] = 0.49910;
812 g[4][1] = -0.28430;
813 //4pi channel
814 g[0][2] = 0;
815 g[1][2] = 0;
816 g[2][2] = 0.70340;
817 g[3][2] = 0.96819;
818 g[4][2] = -0.90100;
819 //eta eta channel
820 g[0][3] = -0.32936;
821 g[1][3] = 0.39910;
822 g[2][3] = 0.22963;
823 g[3][3] = 0.24415;
824 g[4][3] = -0.07252;
825 //eta eta' channel
826 g[0][4] = -0.30906;
827 g[1][4] = 0.31143;
828 g[2][4] = 0.19802;
829 g[3][4] = -0.00522;
830 g[4][4] = 0.17097;
831
832 // Pole masses
833 ma[0] = 0.67460;
834 ma[1] = 1.21094;
835 ma[2] = 1.57896;
836 ma[3] = 1.21900;
837 ma[4] = 1.86602;
838 }
839
840 //Now define the K-matrix pole
841 double rho1sq, rho2sq, rho4sq, rho5sq;
842 EvtComplex rho[5];
843 double f[5][5];
844
845 //Initalize the mass of the resonance
846 const double mpi = 0.13957;
847 const double mK = 0.493677; //using charged K value
848 const double meta = 0.54775; //using PDG value
849 const double metap = 0.95778; //using PDG value
850
851 //Initialize the matrix to value zero
852 EvtComplex K[5][5];
853 for ( int i = 0; i < 5; i++ ) {
854 for ( int j = 0; j < 5; j++ ) {
855 K[i][j] = EvtComplex( 0, 0 );
856 f[i][j] = 0;
857 }
858 }
859
860 //Input the _f[i][j] scattering data
861 double s_scatt = 0.0;
862 if ( solution == 3 )
863 s_scatt = -3.92637;
864 else if ( solution == 1 )
865 s_scatt = -5.0;
866 else if ( solution == 2 )
867 s_scatt = -5.0;
868 const double sa = 1.0;
869 const double sa_0 = -0.15;
870 if ( solution == 3 ) {
871 f[0][0] = 0.23399; // f^scatt
872 f[0][1] = 0.15044;
873 f[0][2] = -0.20545;
874 f[0][3] = 0.32825;
875 f[0][4] = 0.35412;
876 } else if ( solution == 1 ) {
877 f[0][0] = 0.04214; // f^scatt
878 f[0][1] = 0.19865;
879 f[0][2] = -0.63764;
880 f[0][3] = 0.44063;
881 f[0][4] = 0.36717;
882 } else if ( solution == 2 ) {
883 f[0][0] = 0.26447; // f^scatt
884 f[0][1] = 0.10400;
885 f[0][2] = -0.35445;
886 f[0][3] = 0.31596;
887 f[0][4] = 0.42483;
888 }
889 f[1][0] = f[0][1];
890 f[2][0] = f[0][2];
891 f[3][0] = f[0][3];
892 f[4][0] = f[0][4];
893
894 //Now construct the phase-space factor
895 //For eta-eta' there is no difference term
896 rho1sq = 1. - pow( mpi + mpi, 2 ) / s; //pi+ pi- phase factor
897 if ( rho1sq >= 0 )
898 rho[0] = EvtComplex( sqrt( rho1sq ), 0 );
899 else
900 rho[0] = EvtComplex( 0, sqrt( -rho1sq ) );
901
902 rho2sq = 1. - pow( mK + mK, 2 ) / s;
903 if ( rho2sq >= 0 )
904 rho[1] = EvtComplex( sqrt( rho2sq ), 0 );
905 else
906 rho[1] = EvtComplex( 0, sqrt( -rho2sq ) );
907
908 //using the A&S 4pi phase space Factor:
909 //Shit, not continue
910 if ( s <= 1 ) {
911 const double real = 1.2274 + .00370909 / ( s * s ) - .111203 / s -
912 6.39017 * s + 16.8358 * s * s -
913 21.8845 * s * s * s + 11.3153 * s * s * s * s;
914 const double cont32 = sqrt( 1.0 - ( 16.0 * mpi * mpi ) );
915 rho[2] = EvtComplex( cont32 * real, 0 );
916 } else
917 rho[2] = EvtComplex( sqrt( 1. - 16. * mpi * mpi / s ), 0 );
918
919 rho4sq = 1. - pow( meta + meta, 2 ) / s;
920 if ( rho4sq >= 0 )
921 rho[3] = EvtComplex( sqrt( rho4sq ), 0 );
922 else
923 rho[3] = EvtComplex( 0, sqrt( -rho4sq ) );
924
925 rho5sq = 1. - pow( meta + metap, 2 ) / s;
926 if ( rho5sq >= 0 )
927 rho[4] = EvtComplex( sqrt( rho5sq ), 0 );
928 else
929 rho[4] = EvtComplex( 0, sqrt( -rho5sq ) );
930
931 double smallTerm = 1; // Factor to prevent divergences.
932
933 // Check if some pole may arise problems.
934 for ( int pole = 0; pole < 5; pole++ )
935 if ( fabs( pow( ma[pole], 2 ) - s ) < PRECISION )
936 smallTerm = pow( ma[pole], 2 ) - s;
937
938 //now sum all the pole
939 //equation (3) in the E791 K-matrix paper
940 for ( int i = 0; i < 5; i++ ) {
941 for ( int j = 0; j < 5; j++ ) {
942 for ( int pole_index = 0; pole_index < 5; pole_index++ ) {
943 const double A = g[pole_index][i] * g[pole_index][j];
944 const double B = ma[pole_index] * ma[pole_index] - s;
945
946 if ( fabs( B ) < PRECISION )
947 K[i][j] += EvtComplex( A, 0 );
948 else
949 K[i][j] += EvtComplex( A / B, 0 ) * smallTerm;
950 }
951 }
952 }
953
954 //now add the SVT part
955 for ( int i = 0; i < 5; i++ ) {
956 for ( int j = 0; j < 5; j++ ) {
957 const double C = f[i][j] * ( 1.0 - s_scatt );
958 const double D = ( s - s_scatt );
959 K[i][j] += EvtComplex( C / D, 0 ) * smallTerm;
960 }
961 }
962
963 //Fix the bug in the FOCUS paper
964 //Include the Alder zero term:
965 for ( int i = 0; i < 5; i++ ) {
966 for ( int j = 0; j < 5; j++ ) {
967 const double E = ( s - ( sa * mpi * mpi * 0.5 ) ) * ( 1.0 - sa_0 );
968 const double F = ( s - sa_0 );
969 K[i][j] *= EvtComplex( E / F, 0 );
970 }
971 }
972
973 //This is not correct!
974 //(1-ipK) != (1-iKp)
975 static thread_local EvtMatrix<EvtComplex> mat;
976 mat.setRange(
977 5 ); // Try to do in only the first time. DEFINE ALLOCATION IN CONSTRUCTOR.
978
979 for ( int row = 0; row < 5; row++ )
980 for ( int col = 0; col < 5; col++ )
981 mat( row, col ) = ( row == col ) * smallTerm -
982 EvtComplex( 0., 1. ) * K[row][col] * rho[col];
983
984 EvtMatrix<EvtComplex>* matInverse =
985 mat.inverse(); //The 1st row of the inverse matrix. This matrix is {(I-iKp)^-1}_0j
986 vector<EvtComplex> U1j;
987 for ( int j = 0; j < 5; j++ )
988 U1j.push_back( ( *matInverse )[0][j] );
989
990 delete matInverse;
991
992 //this calculates final F0 factor
993 EvtComplex value( 0, 0 );
994 if ( index <= 5 ) {
995 //this calculates the beta_idx Factors
996 for ( int j = 0; j < 5; j++ ) { // sum for 5 channel
997 const EvtComplex top = U1j[j] * g[index - 1][j];
998 const double bottom = ma[index - 1] * ma[index - 1] - s;
999
1000 if ( fabs( bottom ) < PRECISION )
1001 value += top;
1002 else
1003 value += top / bottom * smallTerm;
1004 }
1005 } else {
1006 //this calculates fprod Factors
1007 value += U1j[0];
1008 value += U1j[1] * m_fr12prod;
1009 value += U1j[2] * m_fr13prod;
1010 value += U1j[3] * m_fr14prod;
1011 value += U1j[4] * m_fr15prod;
1012
1013 value *= ( 1 - m_s0prod ) / ( s - m_s0prod ) * smallTerm;
1014 }
1015
1016 return value;
1017}
1018
1019//replace Breit-Wigner with LASS
1021 const EvtTwoBodyVertex& vd ) const
1022{
1023 const double m = kd.m();
1024 const double q = kd.p();
1025 const double GammaM = m_g0 * vd.widthFactor( kd ); // running width;
1026
1027 //calculate the background phase motion
1028 const double cot_deltaB = 1.0 / ( m_a * q ) + 0.5 * m_r * q;
1029 const double deltaB = atan( 1.0 / cot_deltaB );
1030 const double totalB = deltaB + m_phiB;
1031
1032 //calculate the resonant phase motion
1033 const double deltaR = atan( ( m_m0 * GammaM / ( m_m0 * m_m0 - m * m ) ) );
1034 const double totalR = deltaR + m_phiR;
1035
1036 //sum them up
1037 EvtComplex bkgB, resT;
1038 bkgB = EvtComplex( m_Blass * sin( totalB ), 0 ) *
1039 EvtComplex( cos( totalB ), sin( totalB ) );
1040 resT = EvtComplex( m_R * sin( deltaR ), 0 ) *
1041 EvtComplex( cos( totalR ), sin( totalR ) ) *
1042 EvtComplex( cos( 2 * totalB ), sin( 2 * totalB ) );
1043
1044 EvtComplex T;
1045 if ( m_cutoff > 0 && m > m_cutoff )
1046 T = resT;
1047 else
1048 T = bkgB + resT;
1049
1050 if ( m_scaleByMOverQ )
1051 T *= ( m / q );
1052
1053 return T;
1054}
1055
1056EvtComplex EvtDalitzReso::flatte( const double s ) const
1057{
1058 EvtComplex w;
1059
1060 for ( const auto& param : m_flatteParams ) {
1061 const double m1 = param.m1();
1062 const double m2 = param.m2();
1063 const double g = param.g();
1064 w += ( g * g *
1065 sqrtCplx( ( 1 - ( ( m1 - m2 ) * ( m1 - m2 ) ) / ( s ) ) *
1066 ( 1 - ( ( m1 + m2 ) * ( m1 + m2 ) ) / ( s ) ) ) );
1067 }
1068
1069 const EvtComplex denom = m_m0 * m_m0 - s - EvtComplex( 0, 1 ) * w;
1070
1071 return EvtComplex( 1.0, 0.0 ) / denom;
1072}
double real(const EvtComplex &c)
EvtComplex exp(const EvtComplex &c)
#define PRECISION
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=nullptr)
Definition EvtReport.cpp:32
@ EVTGEN_INFO
Definition EvtReport.hh:52
static const double pi
Definition EvtConst.hh:26
static const double twoPi
Definition EvtConst.hh:27
double bigM() const
double m(EvtCyclic3::Index) const
double q(EvtCyclic3::Pair) const
double cosTh(EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes) const
EvtComplex sqrtCplx(const double in) const
EvtComplex propGauss(const double m0, const double s0, const double m) const
std::vector< EvtFlatteParam > m_flatteParams
EvtDalitzPlot m_dp
EvtComplex propBreitWigner(const double m0, const double g0, const double m) const
EvtComplex m_amp_mix
EvtComplex lass(const EvtTwoBodyKine &kd, const EvtTwoBodyVertex &vd) const
EvtComplex propBreitWignerRelCoupled(const double m0, const EvtComplex &g1, const EvtComplex &g2, const double m) const
EvtSpinType::spintype m_spin
EvtComplex propGounarisSakurai(const double m0, const double g0, const double k0, const double m, const double g, const double k) const
EvtComplex numerator(const EvtDalitzPoint &p, const EvtTwoBodyVertex &vb, const EvtTwoBodyVertex &vd, const EvtTwoBodyKine &kb, const EvtTwoBodyKine &kd) const
EvtComplex m_fr13prod
CouplingType m_coupling2
EvtCyclic3::Pair m_pairRes
EvtComplex flatte(const double s) const
double GS_h(const double m, const double k) const
EvtComplex m_fr14prod
double GS_d(const double m0, const double k0) const
double GS_f(const double m0, const double g0, const double k0, const double m, const double k) const
EvtComplex m_fr15prod
EvtComplex mixFactor(const EvtComplex &prop, const EvtComplex &prop_mix) const
double angDep(const EvtDalitzPoint &p) const
EvtComplex propBreitWignerRel(const double m0, const double g0, const double m) const
EvtComplex m_fr12prod
EvtComplex psFactor(const double ma, const double mb, const double m) const
EvtComplex evaluate(const EvtDalitzPoint &p) const
EvtComplex Fvector(const double s, const int index) const
EvtCyclic3::Pair m_pairAng
double GS_dhods(const double m0, const double k0) const
void setRange(int range)
Definition EvtMatrix.hh:52
EvtMatrix * inverse()
Definition EvtMatrix.hh:143
static double getMass(EvtId i)
Definition EvtPDL.cpp:311
static EvtId getId(const std::string &name)
Definition EvtPDL.cpp:283
static int getSpin2(spintype stype)
double m(Index i=AB) const
double p(Index i=AB) const
double formFactor(EvtTwoBodyKine x) const
void set_f(double R)
double pD() const
double widthFactor(EvtTwoBodyKine x) const
double phaseSpaceFactor(EvtTwoBodyKine x, EvtTwoBodyKine::Index) const
static double d(int j, int m1, int m2, double theta)
Index other(Index i, Index j)