BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
D0ToKSpipi.cxx
Go to the documentation of this file.
1#include "D0ToKSpipi.h"
2#include "TMath.h"
3#include <complex>
4#include <iostream>
5#include <math.h>
6#include <stdlib.h>
7#include <string>
8#include <vector>
9
10#include "CLHEP/Matrix/Matrix.h"
11#include "CLHEP/Matrix/SymMatrix.h"
12#include "CLHEP/Matrix/Vector.h"
13#include "CLHEP/Random/RandFlat.h"
14#include "CLHEP/Vector/LorentzVector.h"
15#include "CLHEP/Vector/ThreeVector.h"
16#include "CLHEP/Vector/TwoVector.h"
17using CLHEP::Hep2Vector;
18using CLHEP::Hep3Vector;
19using CLHEP::HepLorentzVector;
20using CLHEP::HepVector;
21using namespace std;
22
23D0ToKSpipi::~D0ToKSpipi() {}
24
25void D0ToKSpipi::init() {
26 // std::cout << "D0ToKSpipi ==> Initialization !" << std::endl;
27
28 _nd = 3;
29 tan2thetaC = ( 0.22650 * 0.22650 ) /
30 ( 1. - ( 0.22650 * 0.22650 ) ); // sin(theta_C) = 0.22650 +/- 0.00048
31 pi180inv = 1.0 * 3.1415926 / 180;
32
33 mass_R[0] = 0.77526;
34 width_R[0] = 0.14740;
35 spin_R[0] = 1;
36 ar[0] = 1;
37 phir[0] = 0;
38 mass_R[1] = 0.78266;
39 width_R[1] = 0.00868;
40 spin_R[1] = 1;
41 ar[1] = 0.037606;
42 phir[1] = 109.677;
43 mass_R[2] = 1.27550;
44 width_R[2] = 0.18670;
45 spin_R[2] = 2;
46 ar[2] = 1.54909;
47 phir[2] = -42.7425;
48 mass_R[3] = 1.46500;
49 width_R[3] = 0.40000;
50 spin_R[3] = 1;
51 ar[3] = 3.70735;
52 phir[3] = 103.644;
53 mass_R[4] = 0.89167;
54 width_R[4] = 0.0514;
55 spin_R[4] = 1;
56 ar[4] = 1.86093;
57 phir[4] = 136.529;
58 mass_R[5] = 1.42730;
59 width_R[5] = 0.10000;
60 spin_R[5] = 2;
61 ar[5] = 1.74288;
62 phir[5] = -48.0968;
63 mass_R[6] = 1.71800;
64 width_R[6] = 0.3220;
65 spin_R[6] = 1;
66 ar[6] = 3.31;
67 phir[6] = -118.2;
68 mass_R[7] = 1.41400;
69 width_R[7] = 0.2320;
70 spin_R[7] = 1;
71 ar[7] = 0.171672;
72 phir[7] = -68.41;
73 mass_R[8] = 0.89167;
74 width_R[8] = 0.0514;
75 spin_R[8] = 1;
76 ar[8] = 0.164;
77 phir[8] = -42.2;
78 mass_R[9] = 1.42730;
79 width_R[9] = 0.1000;
80 spin_R[9] = 2;
81 ar[9] = 0.1;
82 phir[9] = -89.6;
83 mass_R[10] = 1.41400;
84 width_R[10] = 0.2320;
85 spin_R[10] = 1;
86 ar[10] = 0.21;
87 phir[10] = 150.2;
88 mass_R[11] = 1.42500;
89 width_R[11] = 0.2700;
90 spin_R[11] = 1;
91 ar[11] = 2.78276;
92 phir[11] = 97.9608;
93 mass_R[12] = 1.42500;
94 width_R[12] = 0.2700;
95 spin_R[12] = 1;
96 ar[12] = 0.11;
97 phir[12] = 162.3;
98
99 beta[0] = complex<double>( 0.255303 * cos( 47.8861 * pi180inv ),
100 0.255303 * sin( 47.8861 * pi180inv ) );
101 beta[1] = complex<double>( 13.4446 * cos( -5.11127 * pi180inv ),
102 13.4446 * sin( -5.11127 * pi180inv ) );
103 beta[2] = complex<double>( 38.8496 * cos( -30.06 * pi180inv ),
104 38.8496 * sin( -30.06 * pi180inv ) );
105 beta[3] = complex<double>( 13.1086 * cos( -81.4148 * pi180inv ),
106 13.1086 * sin( -81.4148 * pi180inv ) );
107 beta[4] = complex<double>( 0., 0. );
108
109 fprod[0] = complex<double>( 5.08049 * cos( -182.312 * pi180inv ),
110 5.08049 * sin( -182.312 * pi180inv ) );
111 fprod[1] = complex<double>( 17.2388 * cos( -219.209 * pi180inv ),
112 17.2388 * sin( -219.209 * pi180inv ) );
113 fprod[2] = complex<double>( 19.0145 * cos( -76.9884 * pi180inv ),
114 19.0145 * sin( -76.9884 * pi180inv ) );
115 fprod[3] = complex<double>( 11.9875 * cos( -190.502 * pi180inv ),
116 11.9875 * sin( -190.502 * pi180inv ) );
117 fprod[4] = complex<double>( 0., 0. );
118 // beta.push_back( complex<double>( 0.255303*cos( 47.8861 *pi180inv), 0.255303*sin( 47.8861
119 // *pi180inv)) ); beta.push_back( complex<double>(13.4446 *cos( -5.11127*pi180inv), 13.4446
120 // *sin( -5.11127*pi180inv)) ); beta.push_back( complex<double>(38.8496 *cos(-30.06
121 // *pi180inv), 38.8496 *sin(-30.06 *pi180inv)) ); beta.push_back( complex<double>(13.1086
122 // *cos(-81.4148 *pi180inv), 13.1086 *sin(-81.4148 *pi180inv)) ); beta.push_back(
123 // complex<double>(0., 0.) );
124
125 // fprod.push_back(
126 // complex<double>(5.08049*cos(-182.312*pi180inv), 5.08049*sin(-182.312*pi180inv)));
127 // fprod.push_back(
128 // complex<double>(17.2388*cos(-219.209*pi180inv), 17.2388*sin(-219.209*pi180inv)));
129 // fprod.push_back(
130 // complex<double>(19.0145*cos(-76.9884*pi180inv), 19.0145*sin(-76.9884*pi180inv)));
131 // fprod.push_back(
132 // complex<double>(11.9875*cos(-190.502*pi180inv), 11.9875*sin(-190.502*pi180inv)));
133 // fprod.push_back( complex<double>(0., 0.));
134
135 ma[0] = 0.651;
136 g[0][0] = 0.22889;
137 g[0][1] = -0.55377;
138 g[0][2] = 0;
139 g[0][3] = -0.39899;
140 g[0][4] = -0.34639;
141 ma[1] = 1.20360;
142 g[1][0] = 0.94128;
143 g[1][1] = 0.55095;
144 g[1][2] = 0;
145 g[1][3] = 0.39065;
146 g[1][4] = 0.31503;
147 ma[2] = 1.55817;
148 g[2][0] = 0.36856;
149 g[2][1] = 0.23888;
150 g[2][2] = 0.55639;
151 g[2][3] = 0.18340;
152 g[2][4] = 0.18681;
153 ma[3] = 1.21000;
154 g[3][0] = 0.33650;
155 g[3][1] = 0.40907;
156 g[3][2] = 0.85679;
157 g[3][3] = 0.19906;
158 g[3][4] = -0.00984;
159 ma[4] = 1.82206;
160 g[4][0] = 0.18171;
161 g[4][1] = -0.17558;
162 g[4][2] = -0.79658;
163 g[4][3] = -0.00355;
164 g[4][4] = 0.22358;
165
166 // Hadronic parameters for tag modes: 0=no-specified, 1=Kpi, 2=Kpipi0, 3=K3pi
167 rd[0] = 0.0;
168 rd[1] = 0.0586;
169 rd[2] = 0.0440;
170 rd[3] = 0.0546;
171 deltad[0] = 0.0;
172 deltad[1] = 194.7 * pi180inv;
173 deltad[2] = 196.0 * pi180inv;
174 deltad[3] = 167.0 * pi180inv;
175 Rf[0] = 0.0;
176 Rf[1] = 1.0;
177 Rf[2] = 0.78;
178 Rf[3] = 0.52;
179
180 return;
181}
182
183complex<double> D0ToKSpipi::Amp_PFT( vector<double> k0l, vector<double> pip,
184 vector<double> pim ) {
185 // Breit-Wigner lineshapes
186 vector<double> pD;
187 pD.clear();
188 if ( k0l.size() != 4 || pip.size() != 4 || pim.size() != 4 )
189 cout << "ERROR in KSPIPI daughter 4 momentum" << endl;
190 for ( int i = 0; i < k0l.size(); i++ ) { pD.push_back( k0l[i] + pip[i] + pim[i] ); }
191
192 complex<double> DK2piRes0 = Resonance2( pD, pip, pim, ar[0], phir[0], width_R[0], mass_R[0],
193 spin_R[0] ); // ar, phir, width, mass, spin Rho770
194 complex<double> DK2piRes1 = Resonance2( pD, pip, pim, ar[1], phir[1], width_R[1], mass_R[1],
195 spin_R[1] ); // ar, phir, width, mass, spin Omega782
196 complex<double> DK2piRes2 = Resonance2( pD, pip, pim, ar[2], phir[2], width_R[2], mass_R[2],
197 spin_R[2] ); // ar, phir, width, mass, spin ftwo1270
198 complex<double> DK2piRes3 = Resonance2( pD, pip, pim, ar[3], phir[3], width_R[3], mass_R[3],
199 spin_R[3] ); // ar, phir, width, mass, spin Rho1450
200 complex<double> DK2piRes4 = Resonance2( pD, k0l, pim, ar[4], phir[4], width_R[4], mass_R[4],
201 spin_R[4] ); // ar, phir, width, mass, spin Kstar892-
202 complex<double> DK2piRes5 =
203 Resonance2( pD, k0l, pim, ar[5], phir[5], width_R[5], mass_R[5],
204 spin_R[5] ); // ar, phir, width, mass, spin K2star1430-
205 complex<double> DK2piRes6 =
206 Resonance2( pD, k0l, pim, ar[6], phir[6], width_R[6], mass_R[6],
207 spin_R[6] ); // ar, phir, width, mass, spin Kstar1680-
208 complex<double> DK2piRes7 =
209 Resonance2( pD, k0l, pim, ar[7], phir[7], width_R[7], mass_R[7],
210 spin_R[7] ); // ar, phir, width, mass, spin Kstar1410-
211 complex<double> DK2piRes8 = Resonance2( pD, k0l, pip, ar[8], phir[8], width_R[8], mass_R[8],
212 spin_R[8] ); // ar, phir, width, mass, spin Kstar892+
213 complex<double> DK2piRes9 =
214 Resonance2( pD, k0l, pip, ar[9], phir[9], width_R[9], mass_R[9],
215 spin_R[9] ); // ar, phir, width, mass, spin K2star1430+
216 complex<double> DK2piRes10 =
217 Resonance2( pD, k0l, pip, ar[10], phir[10], width_R[10], mass_R[10],
218 spin_R[10] ); // ar, phir, width, mass, spin Kstar1410+
219 // K-matrix for pipi S-wave
220 complex<double> pipi_s_wave = K_matrix( pip, pim );
221 if ( pipi_s_wave == complex<double>( 9999., 9999. ) ) return 1e-20;
222 // LASS parametrization for Kpi S-wave
223 complex<double> kpi_s_wave =
224 amplitude_LASS( k0l, pip, pim, "k0spim", ar[11], phir[11] * pi180inv );
225 complex<double> dcs_kpi_s_wave =
226 amplitude_LASS( k0l, pip, pim, "k0spip", ar[12], phir[12] * pi180inv );
227
228 complex<double> _tmpAmp = DK2piRes0 + DK2piRes1 + DK2piRes2 + DK2piRes3 + pipi_s_wave;
229 // complex<double> TOT_PFT_AMP = DK2piRes0+ DK2piRes1+ DK2piRes2+ DK2piRes3+ DK2piRes4+
230 // DK2piRes5+ DK2piRes6+ DK2piRes7+ DK2piRes8+ DK2piRes9+ DK2piRes10+ pipi_s_wave +
231 // kpi_s_wave+ dcs_kpi_s_wave ;
232 complex<double> TOT_PFT_AMP = _tmpAmp + DK2piRes4 + DK2piRes5 + DK2piRes6 + DK2piRes7 +
233 DK2piRes8 + DK2piRes9 + DK2piRes10 + kpi_s_wave +
234 dcs_kpi_s_wave;
235 // Coherent sum for pure-flavor-tagged amplitudes (PFT)
236 return TOT_PFT_AMP;
237}
238
239complex<double> D0ToKSpipi::Resonance2( vector<double> p4_p, vector<double> p4_d1,
240 vector<double> p4_d2, double mag, double theta,
241 double gamma, double bwm, int spin ) {
242
243 complex<double> ampl;
244
245 // EvtVector4R p4_d3 = p4_p - p4_d1 - p4_d2;
246 HepLorentzVector _p4_p;
247 _p4_p.setX( p4_p[0] );
248 _p4_p.setY( p4_p[1] );
249 _p4_p.setZ( p4_p[2] );
250 _p4_p.setT( p4_p[3] );
251 HepLorentzVector _p4_d1;
252 _p4_d1.setX( p4_d1[0] );
253 _p4_d1.setY( p4_d1[1] );
254 _p4_d1.setZ( p4_d1[2] );
255 _p4_d1.setT( p4_d1[3] );
256 HepLorentzVector _p4_d2;
257 _p4_d2.setX( p4_d2[0] );
258 _p4_d2.setY( p4_d2[1] );
259 _p4_d2.setZ( p4_d2[2] );
260 _p4_d2.setT( p4_d2[3] );
261 HepLorentzVector _p4_d3 = _p4_p - _p4_d1 - _p4_d2;
262
263 double mAB = ( _p4_d1 + _p4_d2 ).invariantMass();
264 double mBC = ( _p4_d2 + _p4_d3 ).invariantMass();
265 double mAC = ( _p4_d1 + _p4_d3 ).invariantMass();
266 double mA = _p4_d1.invariantMass();
267 double mB = _p4_d2.invariantMass();
268 double mD = _p4_p.invariantMass();
269 double mC = _p4_d3.invariantMass();
270
271 double mR = bwm;
272 double gammaR = gamma;
273 double pAB =
274 sqrt( ( ( ( mAB * mAB - mA * mA - mB * mB ) * ( mAB * mAB - mA * mA - mB * mB ) / 4.0 ) -
275 mA * mA * mB * mB ) /
276 ( mAB * mAB ) );
277 double pR =
278 sqrt( ( ( ( mR * mR - mA * mA - mB * mB ) * ( mR * mR - mA * mA - mB * mB ) / 4.0 ) -
279 mA * mA * mB * mB ) /
280 ( mR * mR ) );
281
282 double pD = ( ( ( mD * mD - mR * mR - mC * mC ) * ( mD * mD - mR * mR - mC * mC ) / 4.0 ) -
283 mR * mR * mC * mC ) /
284 ( mD * mD );
285 if ( pD > 0 ) { pD = sqrt( pD ); }
286 else { pD = 0; }
287 double pDAB =
288 sqrt( ( ( ( mD * mD - mAB * mAB - mC * mC ) * ( mD * mD - mAB * mAB - mC * mC ) / 4.0 ) -
289 mAB * mAB * mC * mC ) /
290 ( mD * mD ) );
291 double fR = 1;
292 double fD = 1;
293 int power = 0;
294 switch ( spin )
295 {
296 case 0:
297 fR = 1.0;
298 fD = 1.0;
299 power = 1;
300 break;
301 case 1:
302 fR = sqrt( 1.0 + 1.5 * 1.5 * pR * pR ) / sqrt( 1.0 + 1.5 * 1.5 * pAB * pAB );
303 fD = sqrt( 1.0 + 5.0 * 5.0 * pD * pD ) / sqrt( 1.0 + 5.0 * 5.0 * pDAB * pDAB );
304 power = 3;
305 break;
306 case 2:
307 fR = sqrt( ( 9 + 3 * pow( ( 1.5 * pR ), 2 ) + pow( ( 1.5 * pR ), 4 ) ) /
308 ( 9 + 3 * pow( ( 1.5 * pAB ), 2 ) + pow( ( 1.5 * pAB ), 4 ) ) );
309 fD = sqrt( ( 9 + 3 * pow( ( 5.0 * pD ), 2 ) + pow( ( 5.0 * pD ), 4 ) ) /
310 ( 9 + 3 * pow( ( 5.0 * pDAB ), 2 ) + pow( ( 5.0 * pDAB ), 4 ) ) );
311 power = 5;
312 break;
313 default: cout << "Incorrect spin in D0ToKSpipi::EvtResonance2.cc\n" << endl;
314 }
315
316 double gammaAB = gammaR * pow( pAB / pR, power ) * ( mR / mAB ) * fR * fR;
317 switch ( spin )
318 {
319 case 0:
320 ampl = mag * complex<double>( cos( theta * pi180inv ), sin( theta * pi180inv ) ) * fR *
321 fD / ( mR * mR - mAB * mAB - complex<double>( 0.0, mR * gammaAB ) );
322 break;
323 case 1:
324 ampl = mag * complex<double>( cos( theta * pi180inv ), sin( theta * pi180inv ) ) *
325 ( fR * fD *
326 ( mAC * mAC - mBC * mBC +
327 ( ( mD * mD - mC * mC ) * ( mB * mB - mA * mA ) / ( mAB * mAB ) ) ) /
328 ( mR * mR - mAB * mAB - complex<double>( 0.0, mR * gammaAB ) ) );
329 break;
330 case 2:
331 ampl = mag * complex<double>( cos( theta * pi180inv ), sin( theta * pi180inv ) ) *
332 ( fR * fD / ( mR * mR - mAB * mAB - complex<double>( 0.0, mR * gammaAB ) ) ) *
333 ( pow( ( mBC * mBC - mAC * mAC +
334 ( mD * mD - mC * mC ) * ( mA * mA - mB * mB ) / ( mAB * mAB ) ),
335 2 ) -
336 ( 1.0 / 3.0 ) *
337 ( mAB * mAB - 2 * mD * mD - 2 * mC * mC +
338 pow( ( mD * mD - mC * mC ) / mAB, 2 ) ) *
339 ( mAB * mAB - 2 * mA * mA - 2 * mB * mB +
340 pow( ( mA * mA - mB * mB ) / mAB, 2 ) ) );
341 break;
342 default: cout << "Incorrect spin in D0ToKSpipi::Resonance2.cc\n" << endl;
343 }
344
345 return ampl;
346}
347
348complex<double> D0ToKSpipi::K_matrix( vector<double> p_pip, vector<double> p_pim ) {
349 const double mD0 = 1.86483;
350 const double mKl = 0.49761;
351 const double mPi = 0.13957;
352 bool reject = false;
353
354 HepLorentzVector _p_pip( p_pip[0], p_pip[1], p_pip[2], p_pip[3] );
355 HepLorentzVector _p_pim( p_pim[0], p_pim[1], p_pim[2], p_pim[3] );
356
357 double mAB = ( _p_pip + _p_pim ).m();
358 double s = mAB * mAB;
359
360 complex<double> n11, n12, n13, n14, n15, n21, n22, n23, n24, n25, n31, n32, n33, n34, n35,
361 n41, n42, n43, n44, n45, n51, n52, n53, n54, n55;
362 double rho1sq, rho2sq, rho4sq, rho5sq;
363 complex<double> rho1, rho2, rho3, rho4, rho5;
364 vector<complex<double>> rho;
365 rho.clear();
366 complex<double> pole, SVT, Adler;
367 complex<double> det;
368 // vector< vector< complex<double> > > i;
369 double f[5][5];
370
371 const double mpi = 0.13957;
372 const double mK = 0.493677;
373 const double meta = 0.54775;
374 const double metap = 0.95778;
375
376 // Init matrices and vectors with zeros
377 // vector< vector< complex<double> > > K;
378 complex<double> K[5][5];
379 complex<double> i[5][5];
380 for ( int k = 0; k < 5; k++ )
381 {
382 // vector< complex<double> > _itemp;
383 // vector< complex<double> > _Ktemp;
384 for ( int l = 0; l < 5; l++ )
385 {
386 //_itemp.push_back(complex<double>(0.,0.));
387 //_Ktemp.push_back(complex<double>(0.,0.));
388 i[k][l] = complex<double>( 0., 0. );
389 K[k][l] = complex<double>( 0., 0. );
390 f[k][l] = 0.;
391 }
392 // i.push_back(_itemp);
393 // K.push_back(_Ktemp);
394 // rho.pus_back(0.);
395 }
396
397 // Fill scattering data values
398 double s_scatt = -3.92637;
399 double sa = 1.0;
400 double sa_0 = -0.15;
401
402 // f_scattering (At least one of the two channels must be pi+pi-)
403 f[0][0] = 0.23399;
404 f[0][1] = 0.15044;
405 f[0][2] = -0.20545;
406 f[0][3] = 0.32825;
407 f[0][4] = 0.35412;
408
409 f[1][0] = f[0][1];
410 f[2][0] = f[0][2];
411 f[3][0] = f[0][3];
412 f[4][0] = f[0][4];
413
414 // Compute phase space factors
415 // rho_0
416 rho1sq = ( 1.0 - ( pow( ( mpi + mpi ), 2 ) / s ) );
417 if ( rho1sq >= 0. ) rho1 = complex<double>( sqrt( rho1sq ), 0. );
418 else rho1 = complex<double>( 0., sqrt( -rho1sq ) );
419 rho.push_back( rho1 );
420
421 // rho_1
422 rho2sq = ( 1.0 - ( pow( ( mK + mK ), 2 ) / s ) );
423 if ( rho2sq >= 0. ) rho2 = complex<double>( sqrt( rho2sq ), 0. );
424 else rho2 = complex<double>( 0., sqrt( -rho2sq ) );
425 rho.push_back( rho2 );
426
427 // rho_2
428 rho3 = complex<double>( 0., 0. );
429 if ( s <= 1 )
430 {
431 double real = 1.2274 + 0.00370909 / ( s * s ) - ( 0.111203 ) / (s)-6.39017 * s +
432 16.8358 * s * s - 21.8845 * s * s * s + 11.3153 * s * s * s * s;
433 double cont32 = sqrt( 1.0 - ( 16.0 * mpi * mpi ) );
434 rho3 = complex<double>( cont32 * real, 0. );
435 }
436 else rho3 = complex<double>( sqrt( 1.0 - ( 16.0 * mpi * mpi / s ) ), 0. );
437 rho.push_back( rho3 );
438
439 // rho_3
440 rho4sq = ( 1.0 - ( pow( ( meta + meta ), 2 ) / s ) );
441 if ( rho4sq >= 0. ) rho4 = complex<double>( sqrt( rho4sq ), 0. );
442 else rho4 = complex<double>( 0., sqrt( -rho4sq ) );
443 rho.push_back( rho4 );
444
445 // rho_4
446 rho5sq = ( 1.0 - ( pow( ( meta + metap ), 2 ) / s ) );
447 if ( rho5sq >= 0. ) rho5 = complex<double>( sqrt( rho5sq ), 0. );
448 else rho5 = complex<double>( 0., sqrt( -rho5sq ) );
449 rho.push_back( rho5 );
450
451 // Sum over the poles [Intermediate channel(k) -> pole(pole_index) -> final channel(l)]
452 for ( int k = 0; k < 5; k++ )
453 {
454 for ( int l = 0; l < 5; l++ )
455 {
456 for ( int pole_index = 0; pole_index < 5; pole_index++ )
457 {
458 double A = g[pole_index][k] * g[pole_index][l];
459 double B = ma[pole_index] * ma[pole_index] - s;
460 K[k][l] = K[k][l] + complex<double>( A / B, 0. );
461 }
462 }
463 }
464
465 // Direct scattering term [k -> l]
466 for ( int k = 0; k < 5; k++ )
467 {
468 for ( int l = 0; l < 5; l++ )
469 {
470 double C = f[k][l] * ( 1.0 - s_scatt );
471 double D = ( s - s_scatt );
472 K[k][l] = K[k][l] + complex<double>( C / D, 0. );
473 }
474 }
475
476 // Multiplying the "Adler zero" term
477 for ( int k = 0; k < 5; k++ )
478 {
479 for ( int l = 0; l < 5; l++ )
480 {
481 double E = ( s - ( sa * mpi * mpi * 0.5 ) ) * ( 1.0 - sa_0 );
482 double F = ( s - sa_0 );
483 K[k][l] = K[k][l] * complex<double>( E / F, 0. );
484 }
485 }
486
487 // (1 - i rho K)_ij
488 n11 = complex<double>( 1., 0. ) - complex<double>( 0., 1. ) * K[0][0] * rho[0];
489 n12 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[0][1] * rho[1];
490 n13 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[0][2] * rho[2];
491 n14 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[0][3] * rho[3];
492 n15 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[0][4] * rho[4];
493
494 n21 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[1][0] * rho[0];
495 n22 = complex<double>( 1., 0. ) - complex<double>( 0., 1. ) * K[1][1] * rho[1];
496 n23 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[1][2] * rho[2];
497 n24 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[1][3] * rho[3];
498 n25 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[1][4] * rho[4];
499
500 n31 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[2][0] * rho[0];
501 n32 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[2][1] * rho[1];
502 n33 = complex<double>( 1., 0. ) - complex<double>( 0., 1. ) * K[2][2] * rho[2];
503 n34 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[2][3] * rho[3];
504 n35 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[2][4] * rho[4];
505
506 n41 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[3][0] * rho[0];
507 n42 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[3][1] * rho[1];
508 n43 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[3][2] * rho[2];
509 n44 = complex<double>( 1., 0. ) - complex<double>( 0., 1. ) * K[3][3] * rho[3];
510 n45 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[3][4] * rho[4];
511
512 n51 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[4][0] * rho[0];
513 n52 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[4][1] * rho[1];
514 n53 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[4][2] * rho[2];
515 n54 = complex<double>( 0., 0. ) - complex<double>( 0., 1. ) * K[4][3] * rho[3];
516 n55 = complex<double>( 1., 0. ) - complex<double>( 0., 1. ) * K[4][4] * rho[4];
517
518 // Compute the determinant for inverse [Looks horrible but TMatrixT does not support complex
519 // quantities; python bindings may help, working on it.]
520 det = ( n15 * n24 * n33 * n42 * n51 - n14 * n25 * n33 * n42 * n51 -
521 n15 * n23 * n34 * n42 * n51 + n13 * n25 * n34 * n42 * n51 +
522 n14 * n23 * n35 * n42 * n51 - n13 * n24 * n35 * n42 * n51 -
523 n15 * n24 * n32 * n43 * n51 + n14 * n25 * n32 * n43 * n51 +
524 n15 * n22 * n34 * n43 * n51 - n12 * n25 * n34 * n43 * n51 -
525 n14 * n22 * n35 * n43 * n51 + n12 * n24 * n35 * n43 * n51 +
526 n15 * n23 * n32 * n44 * n51 - n13 * n25 * n32 * n44 * n51 -
527 n15 * n22 * n33 * n44 * n51 + n12 * n25 * n33 * n44 * n51 +
528 n13 * n22 * n35 * n44 * n51 - n12 * n23 * n35 * n44 * n51 -
529 n14 * n23 * n32 * n45 * n51 + n13 * n24 * n32 * n45 * n51 +
530 n14 * n22 * n33 * n45 * n51 - n12 * n24 * n33 * n45 * n51 -
531 n13 * n22 * n34 * n45 * n51 + n12 * n23 * n34 * n45 * n51 -
532 n15 * n24 * n33 * n41 * n52 + n14 * n25 * n33 * n41 * n52 +
533 n15 * n23 * n34 * n41 * n52 - n13 * n25 * n34 * n41 * n52 -
534 n14 * n23 * n35 * n41 * n52 + n13 * n24 * n35 * n41 * n52 +
535 n15 * n24 * n31 * n43 * n52 - n14 * n25 * n31 * n43 * n52 -
536 n15 * n21 * n34 * n43 * n52 + n11 * n25 * n34 * n43 * n52 +
537 n14 * n21 * n35 * n43 * n52 - n11 * n24 * n35 * n43 * n52 -
538 n15 * n23 * n31 * n44 * n52 + n13 * n25 * n31 * n44 * n52 +
539 n15 * n21 * n33 * n44 * n52 - n11 * n25 * n33 * n44 * n52 -
540 n13 * n21 * n35 * n44 * n52 + n11 * n23 * n35 * n44 * n52 +
541 n14 * n23 * n31 * n45 * n52 - n13 * n24 * n31 * n45 * n52 -
542 n14 * n21 * n33 * n45 * n52 + n11 * n24 * n33 * n45 * n52 +
543 n13 * n21 * n34 * n45 * n52 - n11 * n23 * n34 * n45 * n52 +
544 n15 * n24 * n32 * n41 * n53 - n14 * n25 * n32 * n41 * n53 -
545 n15 * n22 * n34 * n41 * n53 + n12 * n25 * n34 * n41 * n53 +
546 n14 * n22 * n35 * n41 * n53 - n12 * n24 * n35 * n41 * n53 -
547 n15 * n24 * n31 * n42 * n53 + n14 * n25 * n31 * n42 * n53 +
548 n15 * n21 * n34 * n42 * n53 - n11 * n25 * n34 * n42 * n53 -
549 n14 * n21 * n35 * n42 * n53 + n11 * n24 * n35 * n42 * n53 +
550 n15 * n22 * n31 * n44 * n53 - n12 * n25 * n31 * n44 * n53 -
551 n15 * n21 * n32 * n44 * n53 + n11 * n25 * n32 * n44 * n53 +
552 n12 * n21 * n35 * n44 * n53 - n11 * n22 * n35 * n44 * n53 -
553 n14 * n22 * n31 * n45 * n53 + n12 * n24 * n31 * n45 * n53 +
554 n14 * n21 * n32 * n45 * n53 - n11 * n24 * n32 * n45 * n53 -
555 n12 * n21 * n34 * n45 * n53 + n11 * n22 * n34 * n45 * n53 -
556 n15 * n23 * n32 * n41 * n54 + n13 * n25 * n32 * n41 * n54 +
557 n15 * n22 * n33 * n41 * n54 - n12 * n25 * n33 * n41 * n54 -
558 n13 * n22 * n35 * n41 * n54 + n12 * n23 * n35 * n41 * n54 +
559 n15 * n23 * n31 * n42 * n54 - n13 * n25 * n31 * n42 * n54 -
560 n15 * n21 * n33 * n42 * n54 + n11 * n25 * n33 * n42 * n54 +
561 n13 * n21 * n35 * n42 * n54 - n11 * n23 * n35 * n42 * n54 -
562 n15 * n22 * n31 * n43 * n54 + n12 * n25 * n31 * n43 * n54 +
563 n15 * n21 * n32 * n43 * n54 - n11 * n25 * n32 * n43 * n54 -
564 n12 * n21 * n35 * n43 * n54 + n11 * n22 * n35 * n43 * n54 +
565 n13 * n22 * n31 * n45 * n54 - n12 * n23 * n31 * n45 * n54 -
566 n13 * n21 * n32 * n45 * n54 + n11 * n23 * n32 * n45 * n54 +
567 n12 * n21 * n33 * n45 * n54 - n11 * n22 * n33 * n45 * n54 +
568 n14 * n23 * n32 * n41 * n55 - n13 * n24 * n32 * n41 * n55 -
569 n14 * n22 * n33 * n41 * n55 + n12 * n24 * n33 * n41 * n55 +
570 n13 * n22 * n34 * n41 * n55 - n12 * n23 * n34 * n41 * n55 -
571 n14 * n23 * n31 * n42 * n55 + n13 * n24 * n31 * n42 * n55 +
572 n14 * n21 * n33 * n42 * n55 - n11 * n24 * n33 * n42 * n55 -
573 n13 * n21 * n34 * n42 * n55 + n11 * n23 * n34 * n42 * n55 +
574 n14 * n22 * n31 * n43 * n55 - n12 * n24 * n31 * n43 * n55 -
575 n14 * n21 * n32 * n43 * n55 + n11 * n24 * n32 * n43 * n55 +
576 n12 * n21 * n34 * n43 * n55 - n11 * n22 * n34 * n43 * n55 -
577 n13 * n22 * n31 * n44 * n55 + n12 * n23 * n31 * n44 * n55 +
578 n13 * n21 * n32 * n44 * n55 - n11 * n23 * n32 * n44 * n55 -
579 n12 * n21 * n33 * n44 * n55 + n11 * n22 * n33 * n44 * n55 );
580
581 if ( det == complex<double>( 0., 0. ) ) reject = true;
582
583 // The 1st row of the inverse matrix [(1-i\rhoK)^-1]_0j
584 i[0][0] = ( n25 * n34 * n43 * n52 - n24 * n35 * n43 * n52 - n25 * n33 * n44 * n52 +
585 n23 * n35 * n44 * n52 + n24 * n33 * n45 * n52 - n23 * n34 * n45 * n52 -
586 n25 * n34 * n42 * n53 + n24 * n35 * n42 * n53 + n25 * n32 * n44 * n53 -
587 n22 * n35 * n44 * n53 - n24 * n32 * n45 * n53 + n22 * n34 * n45 * n53 +
588 n25 * n33 * n42 * n54 - n23 * n35 * n42 * n54 - n25 * n32 * n43 * n54 +
589 n22 * n35 * n43 * n54 + n23 * n32 * n45 * n54 - n22 * n33 * n45 * n54 -
590 n24 * n33 * n42 * n55 + n23 * n34 * n42 * n55 + n24 * n32 * n43 * n55 -
591 n22 * n34 * n43 * n55 - n23 * n32 * n44 * n55 + n22 * n33 * n44 * n55 ) /
592 det;
593
594 i[0][1] = ( -n15 * n34 * n43 * n52 + n14 * n35 * n43 * n52 + n15 * n33 * n44 * n52 -
595 n13 * n35 * n44 * n52 - n14 * n33 * n45 * n52 + n13 * n34 * n45 * n52 +
596 n15 * n34 * n42 * n53 - n14 * n35 * n42 * n53 - n15 * n32 * n44 * n53 +
597 n12 * n35 * n44 * n53 + n14 * n32 * n45 * n53 - n12 * n34 * n45 * n53 -
598 n15 * n33 * n42 * n54 + n13 * n35 * n42 * n54 + n15 * n32 * n43 * n54 -
599 n12 * n35 * n43 * n54 - n13 * n32 * n45 * n54 + n12 * n33 * n45 * n54 +
600 n14 * n33 * n42 * n55 - n13 * n34 * n42 * n55 - n14 * n32 * n43 * n55 +
601 n12 * n34 * n43 * n55 + n13 * n32 * n44 * n55 - n12 * n33 * n44 * n55 ) /
602 det;
603
604 i[0][2] = ( n15 * n24 * n43 * n52 - n14 * n25 * n43 * n52 - n15 * n23 * n44 * n52 +
605 n13 * n25 * n44 * n52 + n14 * n23 * n45 * n52 - n13 * n24 * n45 * n52 -
606 n15 * n24 * n42 * n53 + n14 * n25 * n42 * n53 + n15 * n22 * n44 * n53 -
607 n12 * n25 * n44 * n53 - n14 * n22 * n45 * n53 + n12 * n24 * n45 * n53 +
608 n15 * n23 * n42 * n54 - n13 * n25 * n42 * n54 - n15 * n22 * n43 * n54 +
609 n12 * n25 * n43 * n54 + n13 * n22 * n45 * n54 - n12 * n23 * n45 * n54 -
610 n14 * n23 * n42 * n55 + n13 * n24 * n42 * n55 + n14 * n22 * n43 * n55 -
611 n12 * n24 * n43 * n55 - n13 * n22 * n44 * n55 + n12 * n23 * n44 * n55 ) /
612 det;
613
614 i[0][3] = ( -n15 * n24 * n33 * n52 + n14 * n25 * n33 * n52 + n15 * n23 * n34 * n52 -
615 n13 * n25 * n34 * n52 - n14 * n23 * n35 * n52 + n13 * n24 * n35 * n52 +
616 n15 * n24 * n32 * n53 - n14 * n25 * n32 * n53 - n15 * n22 * n34 * n53 +
617 n12 * n25 * n34 * n53 + n14 * n22 * n35 * n53 - n12 * n24 * n35 * n53 -
618 n15 * n23 * n32 * n54 + n13 * n25 * n32 * n54 + n15 * n22 * n33 * n54 -
619 n12 * n25 * n33 * n54 - n13 * n22 * n35 * n54 + n12 * n23 * n35 * n54 +
620 n14 * n23 * n32 * n55 - n13 * n24 * n32 * n55 - n14 * n22 * n33 * n55 +
621 n12 * n24 * n33 * n55 + n13 * n22 * n34 * n55 - n12 * n23 * n34 * n55 ) /
622 det;
623
624 i[0][4] = ( n15 * n24 * n33 * n42 - n14 * n25 * n33 * n42 - n15 * n23 * n34 * n42 +
625 n13 * n25 * n34 * n42 + n14 * n23 * n35 * n42 - n13 * n24 * n35 * n42 -
626 n15 * n24 * n32 * n43 + n14 * n25 * n32 * n43 + n15 * n22 * n34 * n43 -
627 n12 * n25 * n34 * n43 - n14 * n22 * n35 * n43 + n12 * n24 * n35 * n43 +
628 n15 * n23 * n32 * n44 - n13 * n25 * n32 * n44 - n15 * n22 * n33 * n44 +
629 n12 * n25 * n33 * n44 + n13 * n22 * n35 * n44 - n12 * n23 * n35 * n44 -
630 n14 * n23 * n32 * n45 + n13 * n24 * n32 * n45 + n14 * n22 * n33 * n45 -
631 n12 * n24 * n33 * n45 - n13 * n22 * n34 * n45 + n12 * n23 * n34 * n45 ) /
632 det;
633
634 double s0_prod = -0.07;
635
636 complex<double> value0( 0., 0. );
637 complex<double> value1( 0., 0. );
638
639 // [(1-i\rhoK)^-1]_0j*P_j {P_j: Production vector}
640 for ( int k = 0; k < 5; k++ )
641 {
642 double u1j_re = real( i[0][k] );
643 double u1j_im = imag( i[0][k] );
644 if ( u1j_re == 0. || u1j_im == 0. ) reject = true;
645
646 // Initial state to K-matrix pole couplings * Pole to intermediate channels coupling
647 for ( int pole_index = 0; pole_index < 5; pole_index++ )
648 {
649 complex<double> A = beta[pole_index] * g[pole_index][k];
650 value0 += ( i[0][k] * A ) / ( ma[pole_index] * ma[pole_index] - s );
651 }
652
653 // Direct initial state to intermediate channels couplings
654 value1 += i[0][k] * fprod[k];
655 }
656
657 // Slowly varying polynomial term for the direct coupling
658 value1 *= ( 1. - s0_prod ) / ( s - s0_prod );
659
660 if ( reject == true ) return complex<double>( 9999., 9999. );
661 else return ( value0 + value1 );
662}
663
664complex<double> D0ToKSpipi::amplitude_LASS( vector<double> p_k0l, vector<double> p_pip,
665 vector<double> p_pim, string reso, double A_r,
666 double Phi_r ) {
667 double mR = 1.425;
668 double gammaR = 0.27;
669 double mab2 = 0.0;
670
671 HepLorentzVector _p_k0l( p_k0l[0], p_k0l[1], p_k0l[2], p_k0l[3] );
672 HepLorentzVector _p_pip( p_pip[0], p_pip[1], p_pip[2], p_pip[3] );
673 HepLorentzVector _p_pim( p_pim[0], p_pim[1], p_pim[2], p_pim[3] );
674 if ( reso == "k0spim" ) mab2 = pow( ( _p_k0l + _p_pim ).m(), 2 );
675 else if ( reso == "k0spip" ) mab2 = pow( ( _p_k0l + _p_pip ).m(), 2 );
676 double s = mab2;
677
678 const double mD0 = 1.86483;
679 const double mKl = 0.49761;
680 const double mPi = 0.13957;
681
682 double _a = 0.113;
683 double _r = -33.8;
684 double _R = 1.0;
685 double _F = 0.96;
686 double _phiR = -1.9146;
687 double _phiF = 0.0017;
688 double fR = 1.0; // K*0(1430) has spin zero
689 int power = 1; // Power is 1 for spin zero
690
691 double mAB = sqrt( mab2 );
692 double mA = mKl;
693 double mB = mPi;
694 double mC = mPi;
695 double mD = mD0;
696
697 double pAB =
698 sqrt( ( ( ( mAB * mAB - mA * mA - mB * mB ) * ( mAB * mAB - mA * mA - mB * mB ) / 4.0 ) -
699 mA * mA * mB * mB ) /
700 ( mAB * mAB ) );
701 double q = pAB;
702
703 double pR =
704 sqrt( ( ( ( mR * mR - mA * mA - mB * mB ) * ( mR * mR - mA * mA - mB * mB ) / 4.0 ) -
705 mA * mA * mB * mB ) /
706 ( mR * mR ) );
707 double q0 = pR;
708
709 // Running width.
710 double g = gammaR * pow( q / q0, power ) * ( mR / mAB ) * fR * fR;
711 complex<double> propagator_relativistic_BreitWigner =
712 1. / ( mR * mR - mAB * mAB - complex<double>( 0., mR * g ) );
713
714 // Non-resonant phase shift
715 double cot_deltaF = 1.0 / ( _a * q ) + 0.5 * _r * q;
716 double qcot_deltaF = 1.0 / _a + 0.5 * _r * q * q;
717
718 // Compute resonant part
719 complex<double> expi2deltaF =
720 complex<double>( qcot_deltaF, q ) / complex<double>( qcot_deltaF, -q );
721 complex<double> resonant_term_T =
722 _R * complex<double>( cos( _phiR + 2 * _phiF ), sin( _phiR + 2 * _phiF ) ) *
723 propagator_relativistic_BreitWigner * mR * gammaR * mR / q0 * expi2deltaF;
724
725 // Compute non-resonant part
726 complex<double> non_resonant_term_F = _F * complex<double>( cos( _phiF ), sin( _phiF ) ) *
727 ( cos( _phiF ) + cot_deltaF * sin( _phiF ) ) *
728 sqrt( s ) / complex<double>( qcot_deltaF, -q );
729
730 // Add non-resonant and resonant terms
731 complex<double> LASS_contribution = non_resonant_term_F + resonant_term_T;
732
733 return complex<double>( A_r * cos( Phi_r ), A_r * sin( Phi_r ) ) * LASS_contribution;
734}
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
imag
double imag(const EvtComplex &c)
Definition EvtComplex.hh:207
meta
double meta
Definition EvtPhokharaDef.hh:22
mpi
double mpi
Definition EvtPhokharaDef.hh:17
mPi
double mPi
Definition FarmMonitorAlg.cxx:29
s
XmlRpcServer s
Definition HelloServer.cpp:11
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
q
****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
mD0
const double mD0
Definition MyConst.h:5
C
***************************************************************************************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
D0ToKSpipi::Amp_PFT
complex< double > Amp_PFT(vector< double > k0l, vector< double > pip, vector< double > pim)
Definition D0ToKSpipi.cxx:183
D0ToKSpipi::init
void init()
Definition D0ToKSpipi.cxx:25
D0ToKSpipi::~D0ToKSpipi
virtual ~D0ToKSpipi()
Definition D0ToKSpipi.cxx:23