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