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