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