BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDalitzReso.cc
Go to the documentation of this file.
2/*****************************************************************************
3 * Project: BaBar detector at the SLAC PEP-II B-factory
4 * Package: EvtGenBase
5 * File: $Id: EvtDalitzReso.cc,v 1.1 2009/05/08 01:59:56 pingrg Exp $
6 *
7 * Description:
8 * Class to compute Dalitz amplitudes based on many models that cannot be
9 * handled with EvtResonance.
10 *
11 * Modification history:
12 * Jordi Garra Tic� 2008/07/03 File created
13 *****************************************************************************/
14
15
16#include <assert.h>
17#include <cmath>
18#include <iostream>
19
20#include <stdlib.h>
23#include "EvtGenBase/EvtPDL.hh"
27
30
31#define PRECISION ( 1.e-3 )
32
35
36
37// single Breit-Wigner
38EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
39 EvtSpinType::spintype spin, double m0, double g0, NumType typeN)
40 : _dp(dp),
41 _pairAng(pairAng),
42 _pairRes(pairRes),
43 _spin(spin),
44 _typeN(typeN),
45 _m0(m0),_g0(g0),
46 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
47 _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
48 _g1(-1.),_g2(-1.),_coupling2(Undefined),
49 _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
50 _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
51{
52 _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
53 _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
54 _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
55 _vd.set_f( 1.5 );
56 assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II); // single BW cannot be K-matrix
57}
58
59
60// Breit-Wigner with electromagnetic mass mixing
61EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
62 EvtSpinType::spintype spin, double m0, double g0, NumType typeN,
63 double m0_mix, double g0_mix, double delta_mix, EvtComplex amp_mix)
64 : _dp(dp),
65 _pairAng(pairAng),
66 _pairRes(pairRes),
67 _spin(spin),
68 _typeN(typeN),
69 _m0(m0),_g0(g0),
70 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
71 _m0_mix(m0_mix),_g0_mix(g0_mix),_delta_mix(delta_mix),_amp_mix(amp_mix),
72 _g1(-1.),_g2(-1.),_coupling2(Undefined),
73 _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
74 _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
75{
76 _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
77 _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
78 _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
79 _vd.set_f( 1.5 );
80 // single BW (with electromagnetic mixing) cannot be K-matrix
81 assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II);
82}
83
84// coupled Breit-Wigner
85EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
86 EvtSpinType::spintype spin, double m0, NumType typeN, double g1, double g2, CouplingType coupling2)
87 : _dp(dp),
88 _pairAng(pairAng),
89 _pairRes(pairRes),
90 _spin(spin),
91 _typeN(typeN),
92 _m0(m0),_g0(-1.),
93 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
94 _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
95 _g1(g1),_g2(g2),_coupling2(coupling2),
96 _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
97 _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
98{
99 _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
100 _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
101 _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
102 _vd.set_f( 1.5 );
103 assert(_coupling2 != Undefined);
104 assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II); // coupled BW cannot be K-matrix
105 assert(_typeN != LASS); // coupled BW cannot be LASS
106 assert(_typeN != NBW); // for coupled BW, only relativistic BW
107}
108
109
110// K-Matrix (A&S)
111EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairRes, std::string nameIndex, NumType typeN,
112 EvtComplex fr12prod, EvtComplex fr13prod, EvtComplex fr14prod, EvtComplex fr15prod, double s0prod)
113 : _dp(dp),
114 _pairRes(pairRes),
115 _typeN(typeN),
116 _m0(0.),_g0(0.),
117 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
118 _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
119 _g1(-1.),_g2(-1.),_coupling2(Undefined),
120 _kmatrix_index(-1),_fr12prod(fr12prod),_fr13prod(fr13prod),_fr14prod(fr14prod),_fr15prod(fr15prod),_s0prod(s0prod),
121 _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
122{
123 assert(_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II);
125 if (nameIndex=="Pole1") _kmatrix_index=1;
126 else if (nameIndex=="Pole2") _kmatrix_index=2;
127 else if (nameIndex=="Pole3") _kmatrix_index=3;
128 else if (nameIndex=="Pole4") _kmatrix_index=4;
129 else if (nameIndex=="Pole5") _kmatrix_index=5;
130 else if (nameIndex=="f11prod") _kmatrix_index=6;
131 else assert(0);
132}
133
134
135// LASS parameterization
136EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairRes,
137 double m0, double g0, double a, double r, double B, double phiB, double R, double phiR)
138 : _dp(dp),
139 _pairRes(pairRes),
140 _typeN(LASS),
141 _m0(m0),_g0(g0),
142 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
143 _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
144 _g1(-1.),_g2(-1.),_coupling2(Undefined),
145 _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
146 _a(a),_r(r),_Blass(B),_phiB(phiB),_R(R),_phiR(phiR)
147{
149 _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
150 _vd.set_f( 1.5 ); // Default values for Blatt-Weisskopf factors.
151}
152
153
154
156 : _dp(other._dp),
157 _pairAng(other._pairAng),
158 _pairRes(other._pairRes),
159 _spin(other._spin),
160 _typeN(other._typeN),
161 _m0(other._m0),_g0(other._g0),
162 _vb(other._vb),_vd(other._vd),
163 _massFirst(other._massFirst),_massSecond(other._massSecond),
164 _m0_mix(other._m0_mix),_g0_mix(other._g0_mix),_delta_mix(other._delta_mix),_amp_mix(other._amp_mix),
165 _g1(other._g1),_g2(other._g2),_coupling2(other._coupling2),
166 _kmatrix_index(other._kmatrix_index),
167 _fr12prod(other._fr12prod),_fr13prod(other._fr13prod),_fr14prod(other._fr14prod),_fr15prod(other._fr15prod),
168 _s0prod(other._s0prod),
169 _a(other._a),_r(other._r),_Blass(other._Blass),_phiB(other._phiB),_R(other._R),_phiR(other._phiR)
170{}
171
172
174{}
175
176
178{
179 double m = sqrt(x.q(_pairRes));
180
181 // do use always hash table (speed up fitting)
182 if (_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II)
183 return Fvector( m*m, _kmatrix_index );
184
185 if (_typeN==LASS)
186 return lass(m*m);
187
188 EvtComplex amp(1.0,0.0);
189
190 if (_dp.bigM() != x.bigM()) _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),x.bigM(),_spin);
191 EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
192 EvtTwoBodyKine vd(_massFirst,_massSecond,m);
193
194 EvtComplex prop(0,0);
195 if (_typeN==NBW) {
196 prop = propBreitWigner(_m0,_g0,m);
197 } else if (_typeN==GAUSS_CLEO || _typeN==GAUSS_CLEO_ZEMACH) {
198 prop = propGauss(_m0,_g0,m);
199 } else {
200 if (_coupling2==Undefined) {
201 // single BW
202 double g = (_g0<=0. || _vd.pD()<=0.)? -_g0 : _g0*_vd.widthFactor(vd); // running width
203 if (_typeN==GS_CLEO || _typeN==GS_CLEO_ZEMACH) {
204 // Gounaris-Sakurai (GS)
205 assert(_massFirst==_massSecond);
206 prop = propGounarisSakurai(_m0,fabs(_g0),_vd.pD(),m,g,vd.p());
207 } else {
208 // standard relativistic BW
209 prop = propBreitWignerRel(_m0,g,m);
210 }
211 } else {
212 // coupled width BW
213 EvtComplex G1,G2;
214 switch (_coupling2) {
215 case PicPic: {
216 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
217 static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
218 G2 = _g2*_g2*psFactor(mPic,mPic,m);
219 break;
220 }
221 case PizPiz: {
222 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
223 static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
224 G2 = _g2*_g2*psFactor(mPiz,mPiz,m);
225 break;
226 }
227 case PiPi: {
228 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
229 static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
230 static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
231 G2 = _g2*_g2*psFactor(mPic,mPic,mPiz,mPiz,m);
232 break;
233 }
234 case KcKc: {
235 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
236 static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
237 G2 = _g2*_g2*psFactor(mKc,mKc,m);
238 break;
239 }
240 case KzKz: {
241 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
242 static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
243 G2 = _g2*_g2*psFactor(mKz,mKz,m);
244 break;
245 }
246 case KK: {
247 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
248 static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
249 static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
250 G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
251 break;
252 }
253 case EtaPic: {
254 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
255 static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
256 static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
257 G2 = _g2*_g2*psFactor(mEta,mPic,m);
258 break;
259 }
260 case EtaPiz: {
261 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
262 static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
263 static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
264 G2 = _g2*_g2*psFactor(mEta,mPiz,m);
265 break;
266 }
267 case PicPicKK: {
268 static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
269 //G1 = _g1*_g1*psFactor(mPic,mPic,m);
270 G1 = _g1*psFactor(mPic,mPic,m);
271 static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
272 static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
273 //G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
274 G2 = _g2*psFactor(mKc,mKc,mKz,mKz,m);
275 break;
276 }
277 default:
278 std::cout << "EvtDalitzReso:evaluate(): PANIC, wrong coupling2 state." << std::endl;
279 assert(0);
280 break;
281 }
282 // calculate standard couple BW propagator
283 if (_coupling2 != WA76)
284 prop = _g1*propBreitWignerRelCoupled(_m0,G1,G2,m);
285 }
286 }
287 amp *= prop;
288
289 // Compute form-factors (Blatt-Weisskopf penetration factor)
290 amp *= _vb.formFactor(vb);
291 amp *= _vd.formFactor(vd);
292
293 // Compute numerator (angular distribution)
294 amp *= numerator(x,vb,vd);
295
296 // Compute electromagnetic mass mixing factor
297 if (_m0_mix>0.) {
298 EvtComplex prop_mix;
299 if (_typeN==NBW) {
300 prop_mix = propBreitWigner(_m0_mix,_g0_mix,m);
301 } else {
302 assert(_g1<0.); // running width only
303 double g_mix = _g0_mix*_vd.widthFactor(vd);
304 prop_mix = propBreitWignerRel(_m0_mix,g_mix,m);
305 }
306 amp *= mixFactor(prop,prop_mix);
307 }
308
309 return amp;
310}
311
312
313EvtComplex EvtDalitzReso::psFactor(double & ma, double & mb, double& m)
314{
315 if (m>(ma+mb)) {
316 EvtTwoBodyKine vd(ma,mb,m);
317 return EvtComplex(0,2*vd.p()/m);
318 } else {
319 // analytical continuation
320 double s = m*m;
321 double phaseFactor_analyticalCont = -0.5*(sqrt(4*ma*ma/s-1)+sqrt(4*mb*mb/s-1));
322 return EvtComplex(phaseFactor_analyticalCont,0);
323 }
324}
325
326
327EvtComplex EvtDalitzReso::psFactor(double & ma1,double & mb1, double & ma2, double & mb2, double& m)
328{
329 return 0.5*(psFactor(ma1,mb1,m)+psFactor(ma2,mb2,m));
330}
331
332
333EvtComplex EvtDalitzReso::propGauss(const double& m0, const double& s0, const double& m)
334{
335 // Gaussian
336 double gauss = 1./sqrt(EvtConst::twoPi)/s0*exp(-(m-m0)*(m-m0)/2./(s0*s0));
337 return EvtComplex(gauss,0.);
338}
339
340
341EvtComplex EvtDalitzReso::propBreitWigner(const double& m0, const double& g0, const double& m)
342{
343 // non-relativistic BW
344 return sqrt(g0/EvtConst::twoPi)/(m-m0-EvtComplex(0.0,g0/2.));
345}
346
347
348EvtComplex EvtDalitzReso::propBreitWignerRel(const double& m0, const double& g0, const double& m)
349{
350 // relativistic BW with real width
351 return 1./(m0*m0-m*m-EvtComplex(0.,m0*g0));
352}
353
354
355
356EvtComplex EvtDalitzReso::propBreitWignerRel(const double& m0, const EvtComplex& g0, const double& m)
357{
358 // relativistic BW with complex width
359 return 1./(m0*m0-m*m-EvtComplex(0.,m0)*g0);
360}
361
362
363EvtComplex EvtDalitzReso::propBreitWignerRelCoupled(const double& m0, const EvtComplex& g1, const EvtComplex& g2, const double& m)
364{
365 // relativistic coupled BW
366 return 1./(m0*m0-m*m-(g1+g2));
367}
368
369EvtComplex EvtDalitzReso::propGounarisSakurai(const double& m0, const double& g0, const double& k0,
370 const double& m, const double& g, const double& k)
371{
372 // Gounaris-Sakurai parameterization of pi+pi- P wave. PRD, Vol61, 112002. PRL, Vol21, 244.
373 // 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)
374 return (1.+GS_d(m0,k0)*g0/m0)/(m0*m0-m*m-EvtComplex(0.,m*g)+GS_f(m0,g0,k0,m,k));
375}
376
377
378inline double EvtDalitzReso::GS_f(const double& m0, const double& g0, const double& k0, const double& m, const double& k)
379{
380 // m: sqrt(s)
381 // m0: nominal resonance mass
382 // k: momentum of pion in resonance rest frame (at m)
383 // k0: momentum of pion in resonance rest frame (at nominal resonance mass)
384 return g0*m0*m0/(k0*k0*k0)*( k*k*(GS_h(m,k)-GS_h(m0,k0)) + (m0*m0-m*m)*k0*k0*GS_dhods(m0,k0) );
385}
386
387inline double EvtDalitzReso::GS_h(const double& m, const double& k)
388{return 2./EvtConst::pi*k/m*log((m+2.*k)/(2.*_massFirst)) ;}
389
390inline double EvtDalitzReso::GS_dhods(const double& m0, const double& k0)
391{return GS_h(m0,k0)*( 0.125/(k0*k0) - 0.5/(m0*m0) ) + 0.5/(EvtConst::pi*m0*m0) ;}
392
393inline double EvtDalitzReso::GS_d(const double& m0, const double& k0)
394{return 3./EvtConst::pi*_massFirst*_massFirst/(k0*k0)*log((m0+2.*k0)/(2.*_massFirst)) +
395 m0/(2.*EvtConst::pi*k0) - _massFirst*_massFirst*m0/(EvtConst::pi*k0*k0*k0) ;}
396
397
398EvtComplex EvtDalitzReso::numerator(const EvtDalitzPoint& x, const EvtTwoBodyKine& vb, const EvtTwoBodyKine& vd)
399{
400 EvtComplex ret(0.,0.);
401
402 // Non-relativistic Breit-Wigner
403 if(NBW == _typeN) {
404 ret = angDep(x);
405 }
406
407 // Standard relativistic Zemach propagator
408 else if(RBW_ZEMACH == _typeN) {
409 ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
410 }
411
412 // Standard relativistic Zemach propagator
413 else if(RBW_ZEMACH2 == _typeN) {
415 if(_spin == EvtSpinType::VECTOR) {
416 ret *= -4.;
417 } else if(_spin == EvtSpinType::TENSOR) {
418 ret *= 16./3.;
419 } else if(_spin != EvtSpinType::SCALAR)
420 assert(0);
421 }
422
423 // Kuehn-Santamaria normalization:
424 else if(RBW_KUEHN == _typeN) {
425 ret = _m0*_m0 * angDep(x);
426 }
427
428 // CLEO amplitude
429 else if( ( RBW_CLEO == _typeN ) || ( GS_CLEO == _typeN ) ||
430 ( RBW_CLEO_ZEMACH == _typeN ) || ( GS_CLEO_ZEMACH == _typeN ) ||
431 ( GAUSS_CLEO == _typeN ) || ( GAUSS_CLEO_ZEMACH == _typeN)) {
432
433 Index iA = other(_pairAng); // A = other(BC)
434 Index iB = common(_pairRes,_pairAng); // B = common(AB,BC)
435 Index iC = other(_pairRes); // C = other(AB)
436
437 double M = x.bigM();
438 double mA = x.m(iA);
439 double mB = x.m(iB);
440 double mC = x.m(iC);
441 double qAB = x.q(combine(iA,iB));
442 double qBC = x.q(combine(iB,iC));
443 double qCA = x.q(combine(iC,iA));
444
445 double M2 = M*M;
446 double m02 = ((RBW_CLEO_ZEMACH == _typeN)||(GS_CLEO_ZEMACH == _typeN)||(GAUSS_CLEO_ZEMACH == _typeN))? qAB : _m0*_m0;
447 double mA2 = mA*mA;
448 double mB2 = mB*mB;
449 double mC2 = mC*mC;
450
451 if (_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
452 else if(_spin == EvtSpinType::VECTOR) {
453 ret = qCA - qBC + (M2 - mC2)*(mB2 - mA2)/m02;
454 } else if(_spin == EvtSpinType::TENSOR) {
455 double x1 = qBC - qCA + (M2 - mC2)*(mA2 - mB2)/m02;
456 double x2 = M2 - mC2;
457 double x3 = qAB - 2*M2 - 2*mC2 + x2*x2/m02;
458 double x4 = mA2 - mB2;
459 double x5 = qAB - 2*mB2 - 2*mA2 + x4*x4/m02;
460 ret = x1*x1 - x3*x5/3.;
461 } else assert(0);
462 }
463
464 return ret;
465}
466
467
468double EvtDalitzReso::angDep(const EvtDalitzPoint& x)
469{
470 // Angular dependece for factorizable amplitudes
471 // unphysical cosines indicate we are in big trouble
472 double cosTh = x.cosTh(_pairAng,_pairRes); // angle between common(reso,ang) and other(reso)
473 if(fabs(cosTh) > 1.) {
474 report(INFO,"EvtGen") << "cosTh " << cosTh << std::endl;
475 assert(0);
476 }
477
478 // in units of half-spin
479 return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
480}
481
482
483EvtComplex EvtDalitzReso::mixFactor(EvtComplex prop, EvtComplex prop_mix)
484{
485 double Delta = _delta_mix*(_m0+_m0_mix);
486 return 1/(1-Delta*Delta*prop*prop_mix)*(1+_amp_mix*Delta*prop_mix);
487}
488
489
490
491EvtComplex EvtDalitzReso::Fvector( double s, int index )
492{
493 assert(index>=1 && index<=6);
494
495 //Define the complex coupling constant
496 //The convection is as follow
497 //i=0 --> pi+ pi-
498 //i=1 --> KK
499 //i=2 --> 4pi
500 //i=3 --> eta eta
501 //i=4 --> eta eta'
502 //The first index is the resonace-pole index
503
504 double g[5][5]; // Coupling constants. The first index is the pole index. The second index is the decay channel
505 double ma[5]; // Pole masses. The unit is in GeV
506
507 int solution = (_typeN==K_MATRIX)? 3 : ( (_typeN==K_MATRIX_I)? 1 : ( (_typeN==K_MATRIX_II)? 2 : 0 ) ) ;
508 if (solution==0) { std::cout << "EvtDalitzReso::Fvector() error. Kmatrix solution incorrectly chosen ! " << std::endl; abort(); }
509
510 if (solution == 3 ) {
511
512 // coupling constants
513 //pi+pi- channel
514 g[0][0]=0.22889;
515 g[1][0]=0.94128;
516 g[2][0]=0.36856;
517 g[3][0]=0.33650;
518 g[4][0]=0.18171;
519 //K+K- channel
520 g[0][1]=-0.55377;
521 g[1][1]=0.55095;
522 g[2][1]=0.23888;
523 g[3][1]=0.40907;
524 g[4][1]=-0.17558;
525 //4pi channel
526 g[0][2]=0;
527 g[1][2]=0;
528 g[2][2]=0.55639;
529 g[3][2]=0.85679;
530 g[4][2]=-0.79658;
531 //eta eta channel
532 g[0][3]=-0.39899;
533 g[1][3]=0.39065;
534 g[2][3]=0.18340;
535 g[3][3]=0.19906;
536 g[4][3]=-0.00355;
537 //eta eta' channel
538 g[0][4]=-0.34639;
539 g[1][4]=0.31503;
540 g[2][4]=0.18681;
541 g[3][4]=-0.00984;
542 g[4][4]=0.22358;
543
544 // Pole masses
545 ma[0]=0.651;
546 ma[1]=1.20360;
547 ma[2]=1.55817;
548 ma[3]=1.21000;
549 ma[4]=1.82206;
550
551 } else if (solution == 1) { // solnI.txt
552
553 // coupling constants
554 //pi+pi- channel
555 g[0][0]=0.31896;
556 g[1][0]=0.85963;
557 g[2][0]=0.47993;
558 g[3][0]=0.45121;
559 g[4][0]=0.39391;
560 //K+K- channel
561 g[0][1]=-0.49998;
562 g[1][1]=0.52402;
563 g[2][1]=0.40254;
564 g[3][1]=0.42769;
565 g[4][1]=-0.30860;
566 //4pi channel
567 g[0][2]=0;
568 g[1][2]=0;
569 g[2][2]=1.0;
570 g[3][2]=1.15088;
571 g[4][2]=0.33999;
572 //eta eta channel
573 g[0][3]=-0.21554;
574 g[1][3]=0.38093;
575 g[2][3]=0.21811;
576 g[3][3]=0.22925;
577 g[4][3]=0.06919;
578 //eta eta' channel
579 g[0][4]=-0.18294;
580 g[1][4]=0.23788;
581 g[2][4]=0.05454;
582 g[3][4]=0.06444;
583 g[4][4]=0.32620;
584
585 // Pole masses
586 ma[0]=0.7369;
587 ma[1]=1.24347;
588 ma[2]=1.62681;
589 ma[3]=1.21900;
590 ma[4]=1.74932;
591
592 } else if (solution == 2) { // solnIIa.txt
593
594 // coupling constants
595 //pi+pi- channel
596 g[0][0]=0.26014;
597 g[1][0]=0.95289;
598 g[2][0]=0.46244;
599 g[3][0]=0.41848;
600 g[4][0]=0.01804;
601 //K+K- channel
602 g[0][1]=-0.57849;
603 g[1][1]=0.55887;
604 g[2][1]=0.31712;
605 g[3][1]=0.49910;
606 g[4][1]=-0.28430;
607 //4pi channel
608 g[0][2]=0;
609 g[1][2]=0;
610 g[2][2]=0.70340;
611 g[3][2]=0.96819;
612 g[4][2]=-0.90100;
613 //eta eta channel
614 g[0][3]=-0.32936;
615 g[1][3]=0.39910;
616 g[2][3]=0.22963;
617 g[3][3]=0.24415;
618 g[4][3]=-0.07252;
619 //eta eta' channel
620 g[0][4]=-0.30906;
621 g[1][4]=0.31143;
622 g[2][4]=0.19802;
623 g[3][4]=-0.00522;
624 g[4][4]=0.17097;
625
626 // Pole masses
627 ma[0]=0.67460;
628 ma[1]=1.21094;
629 ma[2]=1.57896;
630 ma[3]=1.21900;
631 ma[4]=1.86602;
632 }
633
634 //Now define the K-matrix pole
635 double rho1sq,rho2sq,rho4sq,rho5sq;
636 EvtComplex rho[5];
637 double f[5][5];
638
639 //Initalize the mass of the resonance
640 double mpi=0.13957;
641 double mK=0.493677; //using charged K value
642 double meta=0.54775; //using PDG value
643 double metap=0.95778; //using PDG value
644
645 //Initialize the matrix to value zero
646 EvtComplex K[5][5];
647 for(int i=0;i<5;i++) {
648 for(int j=0;j<5;j++) {
649 K[i][j]=EvtComplex(0,0);
650 f[i][j]=0;
651 }
652 }
653
654 //Input the _f[i][j] scattering data
655 double s_scatt=0.0 ;
656 if (solution == 3)
657 s_scatt=-3.92637;
658 else if (solution == 1)
659 s_scatt= -5.0 ;
660 else if (solution == 2)
661 s_scatt= -5.0 ;
662 double sa=1.0;
663 double sa_0=-0.15;
664 if (solution == 3) {
665 f[0][0]=0.23399; // f^scatt
666 f[0][1]=0.15044;
667 f[0][2]=-0.20545;
668 f[0][3]=0.32825;
669 f[0][4]=0.35412;
670 }else if (solution == 1) {
671 f[0][0]=0.04214; // f^scatt
672 f[0][1]=0.19865;
673 f[0][2]=-0.63764;
674 f[0][3]=0.44063;
675 f[0][4]=0.36717;
676 }else if (solution == 2) {
677 f[0][0]=0.26447; // f^scatt
678 f[0][1]=0.10400;
679 f[0][2]=-0.35445;
680 f[0][3]=0.31596;
681 f[0][4]=0.42483;
682 }
683 f[1][0]=f[0][1];
684 f[2][0]=f[0][2];
685 f[3][0]=f[0][3];
686 f[4][0]=f[0][4];
687
688 //Now construct the phase-space factor
689 //For eta-eta' there is no difference term
690 rho1sq = 1. - pow( mpi + mpi, 2 ) / s; //pi+ pi- phase factor
691 if( rho1sq >= 0 )
692 rho[ 0 ] = EvtComplex( sqrt( rho1sq ), 0 );
693 else
694 rho[ 0 ] = EvtComplex( 0, sqrt( -rho1sq ) );
695
696 rho2sq = 1. - pow( mK + mK, 2 ) / s;
697 if( rho2sq >= 0 )
698 rho[ 1 ] = EvtComplex( sqrt( rho2sq ), 0 );
699 else
700 rho[ 1 ] = EvtComplex( 0, sqrt( -rho2sq ) );
701
702 //using the A&S 4pi phase space Factor:
703 //Shit, not continue
704 if( s <= 1 )
705 {
706 double real = 1.2274 + .00370909 / ( s * s ) - .111203 / s - 6.39017 * s + 16.8358*s*s - 21.8845*s*s*s + 11.3153*s*s*s*s;
707 double cont32 = sqrt(1.0-(16.0*mpi*mpi));
708 rho[ 2 ] = EvtComplex( cont32 * real, 0 );
709 }
710 else
711 rho[ 2 ] = EvtComplex( sqrt( 1. - 16. * mpi * mpi / s ), 0 );
712
713 rho4sq = 1. - pow( meta + meta, 2 ) / s;
714 if( rho4sq >= 0 )
715 rho[ 3 ] = EvtComplex( sqrt( rho4sq ), 0 );
716 else
717 rho[ 3 ] = EvtComplex( 0, sqrt( -rho4sq ) );
718
719 rho5sq = 1. - pow( meta + metap, 2 ) / s;
720 if( rho5sq >= 0 )
721 rho[ 4 ] = EvtComplex( sqrt( rho5sq ), 0 );
722 else
723 rho[ 4 ] = EvtComplex( 0, sqrt( -rho5sq ) );
724
725 double smallTerm = 1; // Factor to prevent divergences.
726
727 // Check if some pole may arise problems.
728 for ( int pole = 0; pole < 5; pole++ )
729 if ( fabs( pow( ma[ pole ], 2 ) - s ) < PRECISION )
730 smallTerm = pow( ma[ pole ], 2 ) - s;
731
732 //now sum all the pole
733 //equation (3) in the E791 K-matrix paper
734 for(int i=0;i<5;i++) {
735 for(int j=0;j<5;j++) {
736 for (int pole_index=0;pole_index<5;pole_index++) {
737 double A=g[pole_index][i]*g[pole_index][j];
738 double B=ma[pole_index]*ma[pole_index]-s;
739
740 if ( fabs( B ) < PRECISION )
741 K[ i ][ j ] += EvtComplex( A , 0 );
742 else
743 K[ i ][ j ] += EvtComplex( A / B, 0 ) * smallTerm;
744 }
745 }
746 }
747
748 //now add the SVT part
749 for(int i=0;i<5;i++) {
750 for(int j=0;j<5;j++) {
751 double C=f[i][j]*(1.0-s_scatt);
752 double D=(s-s_scatt);
753 K[ i ][ j ] += EvtComplex( C / D, 0 ) * smallTerm;
754 }
755 }
756
757 //Fix the bug in the FOCUS paper
758 //Include the Alder zero term:
759 for(int i=0;i<5;i++) {
760 for(int j=0;j<5;j++) {
761 double E=(s-(sa*mpi*mpi*0.5))*(1.0-sa_0);
762 double F=(s-sa_0);
763 K[ i ][ j ] *= EvtComplex(E/F,0);
764 }
765 }
766
767 //This is not correct!
768 //(1-ipK) != (1-iKp)
769 static EvtMatrix< EvtComplex > mat;
770 mat.setRange( 5 ); // Try to do in only the first time. DEFINE ALLOCATION IN CONSTRUCTOR.
771
772 for ( int row = 0; row < 5; row++ )
773 for ( int col = 0; col < 5; col++ )
774 mat( row, col ) = ( row == col ) * smallTerm - EvtComplex( 0., 1. ) * K[ row ][ col ] * rho[ col ];
775
776
777 EvtMatrix< EvtComplex >* matInverse = mat.inverse(); //The 1st row of the inverse matrix. This matrix is {(I-iKp)^-1}_0j
778 vector< EvtComplex > U1j;
779 for ( int j = 0; j < 5; j++ )
780 U1j.push_back( (*matInverse)[ 0 ][ j ] );
781
782 delete matInverse;
783
784 //this calculates final F0 factor
785 EvtComplex value( 0, 0 );
786 if (index<=5) {
787 //this calculates the beta_idx Factors
788 for(int j=0;j<5;j++) { // sum for 5 channel
789 EvtComplex top = U1j[j]*g[index-1][j];
790 double bottom = ma[index-1]*ma[index-1]-s;
791
792 if ( fabs( bottom ) < PRECISION )
793 value += top;
794 else
795 value += top / bottom * smallTerm;
796 }
797 } else {
798 //this calculates fprod Factors
799 value += U1j[0];
800 value += U1j[1]*_fr12prod;
801 value += U1j[2]*_fr13prod;
802 value += U1j[3]*_fr14prod;
803 value += U1j[4]*_fr15prod;
804
805 value *= (1-_s0prod)/(s-_s0prod) * smallTerm;
806 }
807
808 return value;
809}
810
811
812//replace Breit-Wigner with LASS
813EvtComplex EvtDalitzReso::lass(double s)
814{
815 EvtTwoBodyKine vd(_massFirst,_massSecond, sqrt(s));
816 double q = vd.p();
817 double GammaM = _g0*_vd.widthFactor(vd); // running width;
818
819 //calculate the background phase motion
820 double cot_deltaB = 1.0/(_a*q) + 0.5*_r*q;
821 double deltaB = atan( 1.0/cot_deltaB);
822 double totalB = deltaB + _phiB ;
823
824 //calculate the resonant phase motion
825 double deltaR = atan((_m0*GammaM/(_m0*_m0 - s)));
826 double totalR = deltaR + _phiR ;
827
828 //sum them up
829 EvtComplex bkgB,resT;
830 bkgB = EvtComplex(_Blass*sin(totalB),0)*EvtComplex(cos(totalB),sin(totalB));
831 resT = EvtComplex(_R*sin(deltaR),0)*EvtComplex(cos(totalR),sin(totalR))*EvtComplex(cos(2*totalB),sin(2*totalB));
832 EvtComplex T = bkgB + resT;
833
834 return T;
835}
836
837
838
839
TF1 * g1
Double_t x[10]
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
#define PRECISION
double meta
ostream & report(Severity severity, const char *facility)
Definition: EvtReport.cc:36
@ INFO
Definition: EvtReport.hh:52
const double mpi
Definition: Gam4pikp.cxx:47
XmlRpcServer s
Definition: HelloServer.cpp:11
double sin(const BesAngle a)
double cos(const BesAngle a)
****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
***************************************************************************************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
static const double pi
Definition: EvtConst.hh:28
static const double twoPi
Definition: EvtConst.hh:29
double bigM() const
double m(EvtCyclic3::Index i) const
EvtComplex evaluate(const EvtDalitzPoint &p)
void setRange(int range)
Definition: EvtMatrix.hh:45
EvtMatrix * inverse()
Definition: EvtMatrix.hh:140
static double getMass(EvtId i)
Definition: EvtPDL.hh:46
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
static int getSpin2(spintype stype)
Definition: EvtSpinType.hh:34
double p(Index i=AB) const
double formFactor(EvtTwoBodyKine x) const
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)
Definition: EvtdFunction.cc:30
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
Index common(Pair i, Pair j)
Definition: EvtCyclic3.cc:228
Pair combine(Index i, Index j)
Definition: EvtCyclic3.cc:158
Index second(Pair i)
Definition: EvtCyclic3.cc:206
Index other(Index i, Index j)
Definition: EvtCyclic3.cc:118
Index first(Pair i)
Definition: EvtCyclic3.cc:195
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition: TUtil.h:27