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