BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToPipPipPimPi0.hh
Go to the documentation of this file.
1
2#ifndef EVTDSTOPIPPIPPIMPI0_HH
3#define EVTDSTOPIPPIPPIMPI0_HH
4
5#include "../EvtGenBase/EvtDecayProb.hh"
6
7#include <math.h>
8
9class EvtParticle;
10
11class EvtDsToPipPipPimPi0 : public EvtDecayProb {
12 // must set as pi+ pi+ pi- pi0 in order!
13
14public:
15 EvtDsToPipPipPimPi0() {}
16 virtual ~EvtDsToPipPipPimPi0();
17
18 void getName( std::string& name );
19 EvtDecayBase* clone();
20
21 void init();
22 void initProbMax();
23
24 void decay( EvtParticle* p );
25
26private:
27 static double mass1[11];
28 static double mass2[11];
29 static double width1[11];
30 static double width2[11];
31 static double amp[11];
32 static double phase[11];
33 static double a10gam_arr[3000];
34 static double a1pgam_arr[3000];
35 static double mD;
36 static double metap;
37 static double mKst0;
38 static double GKst0;
39 static double mKstp;
40 static double Gkstp;
41 static double mk0;
42 static double mass_Kaon;
43 static double mass_Ks;
44 static double mass_Eta;
45 static double mass_Pion;
46 static double math_pi;
47 static double mass_Pion2;
48 static double mass_2Pion;
49 static double math_2pi;
50 static double rD2; // 5*5
51 static double rRes1; // 3*3
52 static double rRes2; // 3*3
53 static double g1;
54 static double g2; // K*0(1430)
55 static double GS1;
56 static double GS2;
57 static double GS3; // 1/(2*math_2pi)
58 static double GS4; // mass_Pion2/math_pi
59 static int G[4][4];
60 static int E[4][4][4][4];
61
62 void Com_Multi( double a1[2], double a2[2], double res[2] ) {
63 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
64 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
65 }
66 void Com_Divide( double a1[2], double a2[2], double res[2] ) {
67 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
68 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
69 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
70 }
71 double SCADot( double a1[4], double a2[4] ) {
72 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
73 return _cal;
74 }
75 double barrier( int l, double sa, double sb, double sc, double r, double mass ) {
76 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
77 if ( q < 0 ) q = 1e-16;
78 if ( q < 0 ) q = -q;
79 double z;
80 z = q * r;
81 double sa0;
82 sa0 = mass * mass;
83 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
84 // if(q0 < 0) q0 = 1e-16;
85 if ( q0 < 0 ) q0 = -q0;
86 double z0 = q0 * r;
87 double F = 0.0;
88 if ( l == 0 ) F = 1;
89 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
90 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
91 // if(l == 1) F = sqrt(1./(1+z));
92 // if(l == 2) F = sqrt(1./(9+3*z+z*z));
93 // if(l == 1) F = sqrt((2.0*r)/(1+z));
94 // if(l == 2) F = sqrt((13.0*r*r)/(9+3*z+z*z));
95 return F;
96 }
97 void calt1( double daug1[4], double daug2[4], double t1[4] ) {
98 double p, pq, tmp;
99 double pa[4], qa[4];
100 for ( int i = 0; i < 4; i++ )
101 {
102 pa[i] = daug1[i] + daug2[i];
103 qa[i] = daug1[i] - daug2[i];
104 }
105 p = SCADot( pa, pa );
106 pq = SCADot( pa, qa );
107 tmp = pq / p;
108 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
109 }
110 void calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
111 double p, r;
112 double pa[4], t1[4];
113 calt1( daug1, daug2, t1 );
114 r = SCADot( t1, t1 ) / 3.0;
115 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
116 p = SCADot( pa, pa );
117 for ( int i = 0; i < 4; i++ )
118 {
119 for ( int j = 0; j < 4; j++ )
120 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
121 }
122 }
123 double projector1( double p4[4], int i, int j ) {
124 double result = 0.;
125 result += -G[i][j];
126 result +=
127 p4[i] * p4[j] / ( p4[0] * p4[0] - p4[1] * p4[1] - p4[2] * p4[2] - p4[3] * p4[3] );
128 return result;
129 }
130 double projector2( double p4[4], int i, int j, int k, int l ) {
131 double result = 0.;
132 result += 0.5 * projector1( p4, i, k ) * projector1( p4, j, l );
133 result += 0.5 * projector1( p4, i, l ) * projector1( p4, j, k );
134 result += ( -1. / 3. ) * projector1( p4, i, j ) * projector1( p4, k, l );
135 return result;
136 }
137 double wid( double mass2, double mass, double sa, double sb, double sc, double r2, int l ) {
138 double widm = 0.;
139 double m = sqrt( sa );
140 double tmp = sb - sc;
141 double tmp1 = sa + tmp;
142 double q = 0.25 * tmp1 * tmp1 / sa - sb;
143 if ( q < 0 ) q = 1e-16;
144 if ( q < 0 ) q = -q;
145 double tmp2 = mass2 + tmp;
146 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
147 // if(q0<0) q0 = 1e-16;
148 if ( q0 < 0 ) q0 = -q0;
149 double z = q * r2;
150 double z0 = q0 * r2;
151 double t = q / q0;
152 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
153 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
154 else if ( l == 2 )
155 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
156 return widm;
157 }
158 double widl1( double mass2, double mass, double sa, double sb, double sc, double r2 ) {
159 double widm = 0.;
160 double m = sqrt( sa );
161 double tmp = sb - sc;
162 double tmp1 = sa + tmp;
163 double q = 0.25 * tmp1 * tmp1 / sa - sb;
164 // if(q<0) q = 1e-16;
165 if ( q < 0 ) q = -q;
166 double tmp2 = mass2 + tmp;
167 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
168 // if(q0<0) q0 = 1e-16;
169 if ( q0 < 0 ) q0 = -q0;
170 double z = q * r2;
171 double z0 = q0 * r2;
172 double F = ( 1 + z0 ) / ( 1 + z );
173 double t = q / q0;
174 widm = t * sqrt( t ) * mass / m * F;
175 return widm;
176 }
177 void propagatorRBW( double mass2, double mass, double width, double sa, double sb, double sc,
178 double r2, int l, double prop[2] ) {
179 double a[2], b[2];
180 a[0] = 1;
181 a[1] = 0;
182 b[0] = mass2 - sa;
183 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
184 Com_Divide( a, b, prop );
185 }
186 void propagatorMW( double mass, double sa, int mode, double prop[2] ) {
187 double a[2], b[2];
188 a[0] = 1;
189 a[1] = 0;
190 b[0] = mass * mass - sa;
191 if ( mode == 0 ) { b[1] = -mass * a10gam_arr[int( sqrt( sa ) * 1000 )]; }
192 else if ( mode == 1 ) { b[1] = -mass * a1pgam_arr[int( sqrt( sa ) * 1000 )]; }
193 Com_Divide( a, b, prop );
194 }
195 void Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
196 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
197 if ( q > 0 )
198 {
199 rho[0] = 2 * sqrt( q / sa );
200 rho[1] = 0;
201 }
202 else if ( q < 0 )
203 {
204 rho[0] = 0;
205 rho[1] = 2 * sqrt( -q / sa );
206 }
207 }
208
209 void rhoab( double sa, double sb, double sc, double res[2] ) {
210 double tmp = sa + sb - sc;
211 double q = 0.25 * tmp * tmp / sa - sb;
212 if ( q >= 0 )
213 {
214 res[0] = 2.0 * sqrt( q / sa );
215 res[1] = 0.0;
216 }
217 else
218 {
219 res[0] = 0.0;
220 res[1] = 2.0 * sqrt( -q / sa );
221 }
222 }
223 void propagatorGS( double mass2, double mass, double width, double sa, double sb, double sc,
224 double r2, double prop[2] ) {
225 double a[2], b[2];
226 double tmp = sb - sc;
227 double tmp1 = sa + tmp;
228 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
229 if ( q2 < 0 ) q2 = 1e-16;
230 if ( q2 < 0 ) q2 = -q2;
231
232 double tmp2 = mass2 + tmp;
233 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
234 // if(q02<0) q02 = 1e-16;
235 if ( q02 < 0 ) q02 = -q02;
236
237 double q = sqrt( q2 );
238 double q0 = sqrt( q02 );
239 double m = sqrt( sa );
240 double q03 = q0 * q02;
241 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
242
243 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
244 double h0 = GS1 * q0 / mass * tmp3;
245 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
246 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
247 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
248
249 a[0] = 1.0 + d * width / mass;
250 a[1] = 0.0;
251 b[0] = mass2 - sa + width * f;
252 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
253 Com_Divide( a, b, prop );
254 }
255 void propagator980( double mass, double sx, double* sb, double* sc, double prop[2] ) {
256 double unit[2] = { 1.0 };
257 double ci[2] = { 0, 1 };
258 double rho1[2];
259 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
260 double rho2[2];
261 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
262
263 double gK_f980 = 0.69466, gPi_f980 = 0.165;
264 double tmp1[2] = { gK_f980, 0 };
265 double tmp11[2];
266 double tmp2[2] = { gPi_f980, 0 };
267 double tmp22[2];
268 Com_Multi( tmp1, rho1, tmp11 );
269 Com_Multi( tmp2, rho2, tmp22 );
270 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
271 double tmp31[2];
272 Com_Multi( tmp3, ci, tmp31 );
273 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
274 Com_Divide( unit, tmp4, prop );
275 }
276 double spinf1( double& p1t, double& p1x, double& p1y, double& p1z, double& p2t, double& p2x,
277 double& p2y, double& p2z, double& p3t, double& p3x, double& p3y, double& p3z,
278 double& p4t, double& p4x, double& p4y, double& p4z ) {
279
280 double p1[4] = { p1t, p1x, p1y, p1z };
281 double p2[4] = { p2t, p2x, p2y, p2z };
282 double p3[4] = { p3t, p3x, p3y, p3z };
283 double p4[4] = { p4t, p4x, p4y, p4z };
284 double p14[4] = { p1t + p4t, p1x + p4x, p1y + p4y, p1z + p4z };
285 double p23[4] = { p2t + p3t, p2x + p3x, p2y + p3y, p2z + p3z };
286
287 double t1_14_23[4];
288 double t1_14[4];
289 double t1_23[4];
290
291 calt1( p14, p23, t1_14_23 );
292 calt1( p1, p4, t1_14 );
293 calt1( p2, p3, t1_23 );
294
295 double temp_PDF = 0.;
296 for ( int a = 0; a < 4; a++ )
297 {
298 for ( int b = 0; b < 4; b++ )
299 {
300 if ( a == b ) continue;
301 for ( int c = 0; c < 4; c++ )
302 {
303 if ( a == c ) continue;
304 if ( b == c ) continue;
305 for ( int d = 0; d < 4; d++ )
306 {
307 if ( a == d ) continue;
308 if ( b == d ) continue;
309 if ( c == d ) continue;
310 temp_PDF += G[a][a] * G[b][b] * G[c][c] * G[d][d] * E[a][b][c][d] * t1_14_23[a] *
311 t1_14[b] * t1_23[c] * ( p1[d] + p2[d] + p3[d] + p4[d] );
312 }
313 }
314 }
315 }
316 return temp_PDF;
317 }
318
319 double spinf2( double& p1t, double& p1x, double& p1y, double& p1z, double& p2t, double& p2x,
320 double& p2y, double& p2z, double& p3t, double& p3x, double& p3y, double& p3z,
321 double& p4t, double& p4x, double& p4y, double& p4z ) {
322
323 double p1[4] = { p1t, p1x, p1y, p1z };
324 double p2[4] = { p2t, p2x, p2y, p2z };
325 double p3[4] = { p3t, p3x, p3y, p3z };
326 double p4[4] = { p4t, p4x, p4y, p4z };
327 double p13[4] = { p1t + p3t, p1x + p3x, p1y + p3y, p1z + p3z };
328 double p24[4] = { p2t + p4t, p2x + p4x, p2y + p4y, p2z + p4z };
329
330 double t1_13_24[4];
331 double t1_13[4];
332 double t1_24[4];
333
334 calt1( p13, p24, t1_13_24 );
335 calt1( p1, p3, t1_13 );
336 calt1( p2, p4, t1_24 );
337
338 double temp_PDF = 0.;
339 for ( int a = 0; a < 4; a++ )
340 {
341 for ( int b = 0; b < 4; b++ )
342 {
343 if ( a == b ) continue;
344 for ( int c = 0; c < 4; c++ )
345 {
346 if ( a == c ) continue;
347 if ( b == c ) continue;
348 for ( int d = 0; d < 4; d++ )
349 {
350 if ( a == d ) continue;
351 if ( b == d ) continue;
352 if ( c == d ) continue;
353 temp_PDF += G[a][a] * G[b][b] * G[c][c] * G[d][d] * E[a][b][c][d] * t1_13_24[a] *
354 t1_13[b] * t1_24[c] * ( p1[d] + p2[d] + p3[d] + p4[d] );
355 }
356 }
357 }
358 }
359 return temp_PDF;
360 }
361
362 double spinf3( double& p1t, double& p1x, double& p1y, double& p1z, double& p2t, double& p2x,
363 double& p2y, double& p2z, double& p3t, double& p3x, double& p3y, double& p3z,
364 double& p4t, double& p4x, double& p4y, double& p4z ) {
365
366 double p1[4] = { p1t, p1x, p1y, p1z };
367 double p2[4] = { p2t, p2x, p2y, p2z };
368 double p3[4] = { p3t, p3x, p3y, p3z };
369 double p4[4] = { p4t, p4x, p4y, p4z };
370 double p134[4] = { p1t + p3t + p4t, p1x + p3x + p4x, p1y + p3y + p4y, p1z + p3z + p4z };
371 double p14[4] = { p1t + p4t, p1x + p4x, p1y + p4y, p1z + p4z };
372
373 double s134 =
374 p134[0] * p134[0] - p134[1] * p134[1] - p134[2] * p134[2] - p134[3] * p134[3];
375 double t1_134_2[4];
376 double t1_14_3[4];
377 double t1_14[4];
378
379 calt1( p134, p2, t1_134_2 );
380 calt1( p14, p3, t1_14_3 );
381 calt1( p1, p4, t1_14 );
382
383 double temp_PDF = 0.;
384 for ( int a = 0; a < 4; a++ )
385 {
386 for ( int b = 0; b < 4; b++ )
387 {
388 for ( int c = 0; c < 4; c++ )
389 {
390 if ( b == c ) continue;
391 for ( int d = 0; d < 4; d++ )
392 {
393 if ( b == d ) continue;
394 if ( c == d ) continue;
395 for ( int e = 0; e < 4; e++ )
396 {
397 if ( b == e ) continue;
398 if ( c == e ) continue;
399 if ( d == e ) continue;
400 temp_PDF += G[a][a] * G[b][b] * G[c][c] * G[d][d] * G[e][e] * t1_134_2[a] *
401 ( -G[a][b] + p134[a] * p134[b] / s134 ) * E[b][c][d][e] *
402 t1_14_3[c] * p134[d] * t1_14[e];
403 }
404 }
405 }
406 }
407 }
408 return temp_PDF;
409 }
410
411 double spinf4( double& p1t, double& p1x, double& p1y, double& p1z, double& p2t, double& p2x,
412 double& p2y, double& p2z, double& p3t, double& p3x, double& p3y, double& p3z,
413 double& p4t, double& p4x, double& p4y, double& p4z ) {
414
415 double p1[4] = { p1t, p1x, p1y, p1z };
416 double p2[4] = { p2t, p2x, p2y, p2z };
417 double p3[4] = { p3t, p3x, p3y, p3z };
418 double p4[4] = { p4t, p4x, p4y, p4z };
419 double p234[4] = { p2t + p3t + p4t, p2x + p3x + p4x, p2y + p3y + p4y, p2z + p3z + p4z };
420 double p24[4] = { p2t + p4t, p2x + p4x, p2y + p4y, p2z + p4z };
421
422 double s234 =
423 p234[0] * p234[0] - p234[1] * p234[1] - p234[2] * p234[2] - p234[3] * p234[3];
424 double t1_234_1[4];
425 double t1_24_3[4];
426 double t1_24[4];
427
428 calt1( p234, p1, t1_234_1 );
429 calt1( p24, p3, t1_24_3 );
430 calt1( p2, p4, t1_24 );
431
432 double temp_PDF = 0.;
433 for ( int a = 0; a < 4; a++ )
434 {
435 for ( int b = 0; b < 4; b++ )
436 {
437 for ( int c = 0; c < 4; c++ )
438 {
439 if ( b == c ) continue;
440 for ( int d = 0; d < 4; d++ )
441 {
442 if ( b == d ) continue;
443 if ( c == d ) continue;
444 for ( int e = 0; e < 4; e++ )
445 {
446 if ( b == e ) continue;
447 if ( c == e ) continue;
448 if ( d == e ) continue;
449 temp_PDF += G[a][a] * G[b][b] * G[c][c] * G[d][d] * G[e][e] * t1_234_1[a] *
450 ( -G[a][b] + p234[a] * p234[b] / s234 ) * E[b][c][d][e] *
451 t1_24_3[c] * p234[d] * t1_24[e];
452 }
453 }
454 }
455 }
456 }
457 return temp_PDF;
458 }
459
460 double spinf5( double& p1t, double& p1x, double& p1y, double& p1z, double& p2t, double& p2x,
461 double& p2y, double& p2z, double& p3t, double& p3x, double& p3y, double& p3z,
462 double& p4t, double& p4x, double& p4y, double& p4z ) {
463
464 double p1[4] = { p1t, p1x, p1y, p1z };
465 double p2[4] = { p2t, p2x, p2y, p2z };
466 double p3[4] = { p3t, p3x, p3y, p3z };
467 double p4[4] = { p4t, p4x, p4y, p4z };
468 double p134[4] = { p1t + p3t + p4t, p1x + p3x + p4x, p1y + p3y + p4y, p1z + p3z + p4z };
469 double p13[4] = { p1t + p3t, p1x + p3x, p1y + p3y, p1z + p3z };
470
471 double s134 =
472 p134[0] * p134[0] - p134[1] * p134[1] - p134[2] * p134[2] - p134[3] * p134[3];
473 double t1_134_2[4];
474 double t1_13_4[4];
475 double t1_13[4];
476
477 calt1( p134, p2, t1_134_2 );
478 calt1( p13, p4, t1_13_4 );
479 calt1( p1, p3, t1_13 );
480
481 double temp_PDF = 0.;
482 for ( int a = 0; a < 4; a++ )
483 {
484 for ( int b = 0; b < 4; b++ )
485 {
486 for ( int c = 0; c < 4; c++ )
487 {
488 if ( b == c ) continue;
489 for ( int d = 0; d < 4; d++ )
490 {
491 if ( b == d ) continue;
492 if ( c == d ) continue;
493 for ( int e = 0; e < 4; e++ )
494 {
495 if ( b == e ) continue;
496 if ( c == e ) continue;
497 if ( d == e ) continue;
498 temp_PDF += G[a][a] * G[b][b] * G[c][c] * G[d][d] * G[e][e] * t1_134_2[a] *
499 ( -G[a][b] + p134[a] * p134[b] / s134 ) * E[b][c][d][e] *
500 t1_13_4[c] * p134[d] * t1_13[e];
501 }
502 }
503 }
504 }
505 }
506 return temp_PDF;
507 }
508
509 double spinf6( double& p1t, double& p1x, double& p1y, double& p1z, double& p2t, double& p2x,
510 double& p2y, double& p2z, double& p3t, double& p3x, double& p3y, double& p3z,
511 double& p4t, double& p4x, double& p4y, double& p4z ) {
512
513 double p1[4] = { p1t, p1x, p1y, p1z };
514 double p2[4] = { p2t, p2x, p2y, p2z };
515 double p3[4] = { p3t, p3x, p3y, p3z };
516 double p4[4] = { p4t, p4x, p4y, p4z };
517 double p234[4] = { p2t + p3t + p4t, p2x + p3x + p4x, p2y + p3y + p4y, p2z + p3z + p4z };
518 double p23[4] = { p2t + p3t, p2x + p3x, p2y + p3y, p2z + p3z };
519
520 double s234 =
521 p234[0] * p234[0] - p234[1] * p234[1] - p234[2] * p234[2] - p234[3] * p234[3];
522 double t1_234_1[4];
523 double t1_23_4[4];
524 double t1_23[4];
525
526 calt1( p234, p1, t1_234_1 );
527 calt1( p23, p4, t1_23_4 );
528 calt1( p2, p3, t1_23 );
529
530 double temp_PDF = 0.;
531 for ( int a = 0; a < 4; a++ )
532 {
533 for ( int b = 0; b < 4; b++ )
534 {
535 for ( int c = 0; c < 4; c++ )
536 {
537 if ( b == c ) continue;
538 for ( int d = 0; d < 4; d++ )
539 {
540 if ( b == d ) continue;
541 if ( c == d ) continue;
542 for ( int e = 0; e < 4; e++ )
543 {
544 if ( b == e ) continue;
545 if ( c == e ) continue;
546 if ( d == e ) continue;
547 temp_PDF += G[a][a] * G[b][b] * G[c][c] * G[d][d] * G[e][e] * t1_234_1[a] *
548 ( -G[a][b] + p234[a] * p234[b] / s234 ) * E[b][c][d][e] *
549 t1_23_4[c] * p234[d] * t1_23[e];
550 }
551 }
552 }
553 }
554 }
555 return temp_PDF;
556 }
557
558 double spinf7( double& p1t, double& p1x, double& p1y, double& p1z, double& p2t, double& p2x,
559 double& p2y, double& p2z, double& p3t, double& p3x, double& p3y, double& p3z,
560 double& p4t, double& p4x, double& p4y, double& p4z ) {
561
562 double p1[4] = { p1t, p1x, p1y, p1z };
563 double p2[4] = { p2t, p2x, p2y, p2z };
564 double p3[4] = { p3t, p3x, p3y, p3z };
565 double p4[4] = { p4t, p4x, p4y, p4z };
566 double p134[4] = { p1t + p3t + p4t, p1x + p3x + p4x, p1y + p3y + p4y, p1z + p3z + p4z };
567 double p34[4] = { p3t + p4t, p3x + p4x, p3y + p4y, p3z + p4z };
568
569 double s134 =
570 p134[0] * p134[0] - p134[1] * p134[1] - p134[2] * p134[2] - p134[3] * p134[3];
571 double t1_134_2[4];
572 double t1_34_1[4];
573 double t1_34[4];
574
575 calt1( p134, p2, t1_134_2 );
576 calt1( p34, p1, t1_34_1 );
577 calt1( p3, p4, t1_34 );
578
579 double temp_PDF = 0.;
580 for ( int a = 0; a < 4; a++ )
581 {
582 for ( int b = 0; b < 4; b++ )
583 {
584 for ( int c = 0; c < 4; c++ )
585 {
586 if ( b == c ) continue;
587 for ( int d = 0; d < 4; d++ )
588 {
589 if ( b == d ) continue;
590 if ( c == d ) continue;
591 for ( int e = 0; e < 4; e++ )
592 {
593 if ( b == e ) continue;
594 if ( c == e ) continue;
595 if ( d == e ) continue;
596 temp_PDF += G[a][a] * G[b][b] * G[c][c] * G[d][d] * G[e][e] * t1_134_2[a] *
597 ( -G[a][b] + p134[a] * p134[b] / s134 ) * E[b][c][d][e] *
598 t1_34_1[c] * p134[d] * t1_34[e];
599 }
600 }
601 }
602 }
603 }
604 return temp_PDF;
605 }
606
607 double spinf8( double& p1t, double& p1x, double& p1y, double& p1z, double& p2t, double& p2x,
608 double& p2y, double& p2z, double& p3t, double& p3x, double& p3y, double& p3z,
609 double& p4t, double& p4x, double& p4y, double& p4z ) {
610
611 double p1[4] = { p1t, p1x, p1y, p1z };
612 double p2[4] = { p2t, p2x, p2y, p2z };
613 double p3[4] = { p3t, p3x, p3y, p3z };
614 double p4[4] = { p4t, p4x, p4y, p4z };
615 double p234[4] = { p2t + p3t + p4t, p2x + p3x + p4x, p2y + p3y + p4y, p2z + p3z + p4z };
616 double p34[4] = { p3t + p4t, p3x + p4x, p3y + p4y, p3z + p4z };
617
618 double s234 =
619 p234[0] * p234[0] - p234[1] * p234[1] - p234[2] * p234[2] - p234[3] * p234[3];
620 double t1_234_1[4];
621 double t1_34_2[4];
622 double t1_34[4];
623
624 calt1( p234, p1, t1_234_1 );
625 calt1( p34, p2, t1_34_2 );
626 calt1( p3, p4, t1_34 );
627
628 double temp_PDF = 0.;
629 for ( int a = 0; a < 4; a++ )
630 {
631 for ( int b = 0; b < 4; b++ )
632 {
633 for ( int c = 0; c < 4; c++ )
634 {
635 if ( b == c ) continue;
636 for ( int d = 0; d < 4; d++ )
637 {
638 if ( b == d ) continue;
639 if ( c == d ) continue;
640 for ( int e = 0; e < 4; e++ )
641 {
642 if ( b == e ) continue;
643 if ( c == e ) continue;
644 if ( d == e ) continue;
645 temp_PDF += G[a][a] * G[b][b] * G[c][c] * G[d][d] * G[e][e] * t1_234_1[a] *
646 ( -G[a][b] + p234[a] * p234[b] / s234 ) * E[b][c][d][e] *
647 t1_34_2[c] * p234[d] * t1_34[e];
648 }
649 }
650 }
651 }
652 }
653 return temp_PDF;
654 }
655};
656
657#endif
p2
double p2[4]
Definition BabayagaNLO.cxx:41
p1
double p1[4]
Definition BabayagaNLO.cxx:41
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
g1
TF1 * g1
Definition DataBase/tau_mode.c:58
mass_Pion
const double mass_Pion
Definition EvtD0ToKSpipipi0.cc:1561
unit
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
Definition FoamA.h:90
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
EvtDsToPipPipPimPi0::getName
void getName(std::string &name)
Definition EvtDsToPipPipPimPi0.cc:31
EvtDsToPipPipPimPi0::decay
void decay(EvtParticle *p)
Definition EvtDsToPipPipPimPi0.cc:161
mk0
const double mk0
Definition inclks.cxx:33
t
int t()
Definition t.c:1