BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDTopipi0pi0.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtDTopipi0pi0.cc
11// the necessary file: EvtDTopipi0pi0.hh
12//
13// Description: D+ -> pi+ pi0 pi0 (BAM-00744)
14//
15// Modification history:
16//
17// Liaoyuan Dong Sat Dec 31 13:54:49 2022 Module created
18// Tue Jan 2 12:12:00 2024 Module updated
19//------------------------------------------------------------------------
20#include "EvtDTopipi0pi0.hh"
30#include "TComplex.h"
31#include "TMatrixD.h"
32#include <cmath>
33#include <fstream>
34#include <stdlib.h>
35#include <string>
36using std::endl;
37
39
40void EvtDTopipi0pi0::getName( std::string& model_name ) { model_name = "DTopipi0pi0"; }
41
43
45 // check that there are 0 arguments
46 checkNArg( 0 );
47 checkNDaug( 3 );
52
53 phi[0] = 9.32300;
54 phi[1] = 0.0;
55 phi[2] = 2.61960;
56 phi[3] = 0.0;
57
58 rho[0] = 1.26430;
59 rho[1] = 1.0;
60 rho[2] = 0.95403;
61 rho[3] = 1.0;
62
63 modetype[0] = 23;
64 modetype[1] = 12;
65 modetype[2] = 12;
66 modetype[3] = 23;
67
68 /*
69 std::cout << "Initializing EvtDTopipi0pi0" << std::endl;
70 for (int i=0; i<4; i++) {
71 cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
72 }
73 */
74 width[0] = 0.1867;
75 width[1] = 0.1502;
76 width[2] = 0.4;
77 width[3] = 0.020;
78
79 mass[0] = 1.2755;
80 mass[1] = 0.7665;
81 mass[2] = 1.465;
82 mass[3] = 0.900;
83
84 mD = 1.86486;
85 mDs = 1.9683;
86 rRes = 9.0;
87 rD = 5.0;
88 metap = 0.95778;
89 mkstr = 0.89594;
90 mk0 = 0.497614;
91 mass_Kaon = 0.49368;
92 mass_Pion = 0.13957;
93 mass_Pi0 = 0.1349766;
94 math_pi = 3.1415926;
95 ma0 = 0.99;
96 Ga0 = 0.0756;
97 meta = 0.547862;
98 mass_Kaon = 0.49368;
99 mass_Eta = 0.547862;
100 mass_Etap = 0.95778;
101
102 GS1 = 0.636619783;
103 GS2 = 0.01860182466;
104 GS3 = 0.1591549458; // 1/(2*math_2pi)
105 GS4 = 0.00620060822; // mass_Pion2/math_pi
106
107 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
108 for ( int i = 0; i < 4; i++ )
109 {
110 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
111 }
112}
113
115
117
119 EvtVector4R D1 = p->getDaug( 0 )->getP4();
120 EvtVector4R D2 = p->getDaug( 1 )->getP4();
121 EvtVector4R D3 = p->getDaug( 2 )->getP4();
122
123 double P1[4], P2[4], P3[4];
124 P1[0] = D1.get( 0 );
125 P1[1] = D1.get( 1 );
126 P1[2] = D1.get( 2 );
127 P1[3] = D1.get( 3 );
128 P2[0] = D2.get( 0 );
129 P2[1] = D2.get( 1 );
130 P2[2] = D2.get( 2 );
131 P2[3] = D2.get( 3 );
132 P3[0] = D3.get( 0 );
133 P3[1] = D3.get( 1 );
134 P3[2] = D3.get( 2 );
135 P3[3] = D3.get( 3 );
136
137 double value;
138 int g0[4] = { 1, 2, 2, 6 };
139 int spin[4] = { 2, 1, 1, 0 };
140 int nstates = 4;
141 double r0[4] = { 3, 3, 3, 3 };
142 double r1[4] = { 5, 5, 5, 5 };
143
144 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value, r0, r1 );
145
146 setProb( value );
147 return;
148}
149
150void EvtDTopipi0pi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
151 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
152 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
153}
154void EvtDTopipi0pi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
155 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
156 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
157 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
158}
159
160double EvtDTopipi0pi0::CalRho4pi( double_t s ) {
161 if ( s >= 1. ) { return sqrt( ( s - 16. * mass_Pion * mass_Pion ) / s ); }
162 else
163 {
164 double_t s0 = 1.2274 + 0.00370909 / ( s * s ) - ( 0.111203 ) / (s)-6.39017 * s +
165 16.8358 * s * s - 21.8845 * s * s * s + 11.3153 * s * s * s * s;
166 double_t gam = s0 * sqrt( 1.0 - ( 16.0 * mass_Pion * mass_Pion ) );
167
168 return gam;
169 }
170}
171
172//------------base---------------------------------
173double EvtDTopipi0pi0::SCADot( double a1[4], double a2[4] ) {
174 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
175 return _cal;
176}
177double EvtDTopipi0pi0::barrier( int l, double sa, double sb, double sc, double r,
178 double mass ) {
179 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
180 if ( q < 0 ) q = 1e-16;
181 double z;
182 z = q * r * r;
183 double sa0;
184 sa0 = mass * mass;
185 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
186 if ( q0 < 0 ) q0 = 1e-16;
187 double z0 = q0 * r * r;
188 double F = 0.0;
189 if ( l == 0 ) F = 1;
190 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
191 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
192 return F;
193}
194//------------------spin-------------------------------------------
195void EvtDTopipi0pi0::calt1( double daug1[4], double daug2[4], double t1[4] ) {
196 double p, pq, tmp;
197 double pa[4], qa[4];
198 for ( int i = 0; i < 4; i++ )
199 {
200 pa[i] = daug1[i] + daug2[i];
201 qa[i] = daug1[i] - daug2[i];
202 }
203 p = SCADot( pa, pa );
204 pq = SCADot( pa, qa );
205 tmp = pq / p;
206 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
207}
208void EvtDTopipi0pi0::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
209 double p, r;
210 double pa[4], t1[4];
211 calt1( daug1, daug2, t1 );
212 r = SCADot( t1, t1 ) / 3.0;
213 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
214 p = SCADot( pa, pa );
215 for ( int i = 0; i < 4; i++ )
216 {
217 for ( int j = 0; j < 4; j++ )
218 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
219 }
220}
221//-------------------prop--------------------------------------------
222void EvtDTopipi0pi0::propagator( double mass2, double mass, double width, double sx,
223 double prop[2] ) {
224 double a[2], b[2];
225 a[0] = 1;
226 a[1] = 0;
227 b[0] = mass2 - sx;
228 b[1] = -mass * width;
229 Com_Divide( a, b, prop );
230}
231double EvtDTopipi0pi0::wid( double mass2, double mass, double sa, double sb, double sc,
232 double r2, int l ) {
233 double widm = 0.;
234 double m = sqrt( sa );
235 double tmp = sb - sc;
236 double tmp1 = sa + tmp;
237 double q = 0.25 * tmp1 * tmp1 / sa - sb;
238 if ( q < 0 ) q = 1e-16;
239 double tmp2 = mass2 + tmp;
240 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
241 if ( q0 < 0 ) q0 = 1e-16;
242 double z = q * r2;
243 double z0 = q0 * r2;
244 double t = q / q0;
245 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
246 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
247 else if ( l == 2 )
248 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
249 return widm;
250}
251double EvtDTopipi0pi0::widl1( double mass2, double mass, double sa, double sb, double sc,
252 double r2 ) {
253 double widm = 0.;
254 double m = sqrt( sa );
255 double tmp = sb - sc;
256 double tmp1 = sa + tmp;
257 double q = 0.25 * tmp1 * tmp1 / sa - sb;
258 if ( q < 0 ) q = 1e-16;
259 double tmp2 = mass2 + tmp;
260 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
261 if ( q0 < 0 ) q0 = 1e-16;
262 double z = q * r2;
263 double z0 = q0 * r2;
264 double F = ( 1 + z0 ) / ( 1 + z );
265 double t = q / q0;
266 widm = t * sqrt( t ) * mass / m * F;
267 return widm;
268}
269void EvtDTopipi0pi0::propagatorRBW( double mass2, double mass, double width, double sa,
270 double sb, double sc, double r2, int l, double prop[2] ) {
271
272 double a[2], b[2];
273 a[0] = 1;
274 a[1] = 0;
275 b[0] = mass2 - sa;
276 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
277 Com_Divide( a, b, prop );
278}
279
280void EvtDTopipi0pi0::propagatorFlatte( double mass, double width, double sa, double sb,
281 double sc, int r, double prop[2] ) {
282 double q, qKsK;
283 double rhoab[2], rhoKsK[2];
284 q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
285 if ( r == 0 ) qKsK = 0.25 * sa - 0.49368 * 0.49368;
286 if ( r == 1 )
287 qKsK =
288 0.25 * ( sa + 3.899750596e-03 ) * ( sa + 3.899750596e-03 ) / sa - 0.497614 * 0.497614;
289 if ( q > 0 )
290 {
291 rhoab[0] = 2 * sqrt( q / sa );
292 rhoab[1] = 0;
293 }
294 if ( q < 0 )
295 {
296 rhoab[0] = 0;
297 rhoab[1] = 2 * sqrt( -q / sa );
298 }
299 if ( qKsK > 0 )
300 {
301 rhoKsK[0] = 2 * sqrt( qKsK / sa );
302 rhoKsK[1] = 0;
303 }
304 if ( qKsK < 0 )
305 {
306 rhoKsK[0] = 0;
307 rhoKsK[1] = 2 * sqrt( -qKsK / sa );
308 }
309 double a[2], b[2];
310 a[0] = 1;
311 a[1] = 0;
312 b[0] = mass * mass - sa + 0.341 * rhoab[1] + 0.892 * 0.341 * rhoKsK[1];
313 b[1] = -( 0.341 * rhoab[0] + 0.892 * 0.341 * rhoKsK[0] );
314 Com_Divide( a, b, prop );
315}
316
317void EvtDTopipi0pi0::rhoab( double sa, double sb, double sc, double res[2] ) {
318 double tmp = sa + sb - sc;
319 double q = 0.25 * tmp * tmp / sa - sb;
320 if ( q >= 0 )
321 {
322 res[0] = 2.0 * sqrt( q / sa );
323 res[1] = 0.0;
324 }
325 else
326 {
327 res[0] = 0.0;
328 res[1] = 2.0 * sqrt( -q / sa );
329 }
330}
331void EvtDTopipi0pi0::rho4Pi( double sa, double res[2] ) {
332 double temp = 1.0 - 0.3116765584 / sa; // 0.3116765584=0.13957*0.13957*16
333 if ( temp >= 0 )
334 {
335 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
336 res[1] = 0.0;
337 }
338 else
339 {
340 res[0] = 0.0;
341 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) ); //?????????????????????
342 }
343}
344
345void EvtDTopipi0pi0::propagatorsigma500( double sa, double sb, double sc, double prop[2] ) {
346 double f = 0.5843 + 1.6663 * sa;
347 const double M = 0.9264; // M(f0500)
348 const double mass2 = 0.85821696; // M*M
349 const double mpi2d2 = 0.00973989245; // mass_Pion2/2 = 0.0194797849/2
350 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
351 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
352 rhoab( sa, sb, sc, rho1s );
353 rhoab( mass2, sb, sc, rho1M );
354 rho4Pi( sa, rho2s );
355 rho4Pi( mass2, rho2M );
356 Com_Divide( rho1s, rho1M, rho1 );
357 Com_Divide( rho2s, rho2M, rho2 );
358 double a[2], b[2];
359 a[0] = 1.0;
360 a[1] = 0.0;
361 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
362 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
363 Com_Divide( a, b, prop );
364}
365void EvtDTopipi0pi0::Flatte_rhoab( double sa, double sb, double rho[2] ) {
366 double q = 1.0 - ( 4 * sb / sa );
367
368 if ( q > 0 )
369 {
370 rho[0] = sqrt( q );
371 rho[1] = 0;
372 }
373 else if ( q < 0 )
374 {
375 rho[0] = 0;
376 rho[1] = sqrt( -q );
377 }
378}
379
380void EvtDTopipi0pi0::propagator980( double mass, double sx, double* sb, double prop[2] ) {
381
382 double gpipi1 = 2.0 / 3.0 * 0.165;
383 double gpipi2 = 1.0 / 3.0 * 0.165;
384 double gKK1 = 0.5 * 0.69465;
385 double gKK2 = 0.5 * 0.69465;
386
387 double unit[2] = { 1.0 };
388 double ci[2] = { 0, 1 };
389 double rho1[2];
390 Flatte_rhoab( sx, sb[0], rho1 );
391 double rho2[2];
392 Flatte_rhoab( sx, sb[1], rho2 );
393 double rho3[2];
394 Flatte_rhoab( sx, sb[2], rho3 );
395 double rho4[2];
396 Flatte_rhoab( sx, sb[3], rho4 );
397
398 double tmp1[2] = { gpipi1, 0 };
399 double tmp11[2];
400 double tmp2[2] = { gpipi2, 0 };
401 double tmp22[2];
402 double tmp3[2] = { gKK1, 0 };
403 double tmp33[2];
404 double tmp4[2] = { gKK2, 0 };
405 double tmp44[2];
406
407 Com_Multi( tmp1, rho1, tmp11 );
408 Com_Multi( tmp2, rho2, tmp22 );
409 Com_Multi( tmp3, rho3, tmp33 );
410 Com_Multi( tmp4, rho4, tmp44 );
411
412 double tmp5[2] = { tmp11[0] + tmp22[0] + tmp33[0] + tmp44[0],
413 tmp11[1] + tmp22[1] + tmp33[1] + tmp44[1] };
414 double tmp51[2];
415 Com_Multi( tmp5, ci, tmp51 );
416 double tmp6[2] = { mass * mass - sx - tmp51[0], -1.0 * tmp51[1] };
417
418 Com_Divide( unit, tmp6, prop );
419}
420
421//------------GS---used by rho----------------------------
422void EvtDTopipi0pi0::propagatorGS( double mass2, double mass, double width, double sa,
423 double sb, double sc, double r2, double prop[2] ) {
424 double a[2], b[2];
425 double tmp = sb - sc;
426 double tmp1 = sa + tmp;
427 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
428 if ( q2 < 0 ) q2 = 1e-16;
429
430 double tmp2 = mass2 + tmp;
431 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
432 if ( q02 < 0 ) q02 = 1e-16;
433
434 double q = sqrt( q2 );
435 double q0 = sqrt( q02 );
436 double m = sqrt( sa );
437 double q03 = q0 * q02;
438 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
439
440 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
441 double h0 = GS1 * q0 / mass * tmp3;
442 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
443 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
444 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
445
446 a[0] = 1.0 + d * width / mass;
447 a[1] = 0.0;
448 b[0] = mass2 - sa + width * f;
449 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
450 Com_Divide( a, b, prop );
451}
452
453double EvtDTopipi0pi0::DDalitz( double P1[4], double P2[4], double P3[4], int Ang, double mass,
454 double rRES0, double rRES1 ) {
455 double pR[4], pD[4];
456 double temp_PDF, v_re;
457 temp_PDF = 0.0;
458 v_re = 0.0;
459 double B[2], s1, s2, s3, sR, sD;
460 for ( int i = 0; i < 4; i++ )
461 {
462 pR[i] = P1[i] + P2[i];
463 pD[i] = pR[i] + P3[i];
464 }
465 s1 = SCADot( P1, P1 );
466 s2 = SCADot( P2, P2 );
467 s3 = SCADot( P3, P3 );
468 sR = SCADot( pR, pR );
469 sD = SCADot( pD, pD );
470 int G[4][4];
471 for ( int i = 0; i != 4; i++ )
472 {
473 for ( int j = 0; j != 4; j++ )
474 {
475 if ( i == j )
476 {
477 if ( i == 0 ) G[i][j] = 1;
478 else G[i][j] = -1;
479 }
480 else G[i][j] = 0;
481 }
482 }
483 if ( Ang == 0 )
484 {
485 B[0] = 1;
486 B[1] = 1;
487 temp_PDF = 1;
488 }
489 if ( Ang == 1 )
490 {
491 B[0] = barrier( 1, sR, s1, s2, rRES0, mass );
492 B[1] = barrier( 1, sD, sR, s3, rRES1, 1.86966 );
493 double t1[4], T1[4];
494 calt1( P1, P2, t1 );
495 calt1( pR, P3, T1 );
496 temp_PDF = 0;
497 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
498 }
499 if ( Ang == 2 )
500 {
501 B[0] = barrier( 2, sR, s1, s2, rRES0, mass );
502 B[1] = barrier( 2, sD, sR, s3, rRES1, 1.86966 );
503 double t2[4][4], T2[4][4];
504 calt2( P1, P2, t2 );
505 calt2( pR, P3, T2 );
506 temp_PDF = 0;
507 for ( int i = 0; i < 4; i++ )
508 {
509 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
510 }
511 }
512 v_re = temp_PDF * B[0] * B[1];
513 return v_re;
514}
515
516TComplex EvtDTopipi0pi0::ResonanceSkm( Double_t& m2 ) {
517 Double_t g11 = 0.22889, g12 = -0.55377, g13 = 0.00, g14 = -0.39899, g15 = -0.34639;
518 Double_t g21 = 0.94128, g22 = 0.55095, g23 = 0.00, g24 = 0.39065, g25 = 0.31503;
519 Double_t g31 = 0.36856, g32 = 0.23888, g33 = 0.55639, g34 = 0.18340, g35 = 0.18681;
520 Double_t g41 = 0.33650, g42 = 0.40907, g43 = 0.85679, g44 = 0.19906, g45 = -0.00984;
521 Double_t g51 = 0.18171, g52 = -0.17558, g53 = -0.79658, g54 = -0.00355, g55 = 0.22358;
522
523 Double_t fs11 = 0.23399, fs12 = 0.15044, fs13 = -0.20545, fs14 = 0.32825, fs15 = 0.35412;
524 Double_t m_1 = 0.65100, m_2 = 1.20360, m_3 = 1.55817, m_4 = 1.21000, m_5 = 1.82206;
525 Double_t ss0 = -3.92637;
526 Double_t sp0 = -0.07; // v1
527 double_t sA0 = -0.15;
528 double_t sA = 1.0;
529
530 double_t km11 = ( g11 * g11 / ( m_1 * m_1 - m2 ) + g21 * g21 / ( m_2 * m_2 - m2 ) +
531 g31 * g31 / ( m_3 * m_3 - m2 ) + g41 * g41 / ( m_4 * m_4 - m2 ) +
532 g51 * g51 / ( m_5 * m_5 - m2 ) + fs11 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
533 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
534 double_t km12 = ( g11 * g12 / ( m_1 * m_1 - m2 ) + g21 * g22 / ( m_2 * m_2 - m2 ) +
535 g31 * g32 / ( m_3 * m_3 - m2 ) + g41 * g42 / ( m_4 * m_4 - m2 ) +
536 g51 * g52 / ( m_5 * m_5 - m2 ) + fs12 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
537 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
538 double_t km13 = ( g11 * g13 / ( m_1 * m_1 - m2 ) + g21 * g23 / ( m_2 * m_2 - m2 ) +
539 g31 * g33 / ( m_3 * m_3 - m2 ) + g41 * g43 / ( m_4 * m_4 - m2 ) +
540 g51 * g53 / ( m_5 * m_5 - m2 ) + fs13 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
541 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
542 double_t km14 = ( g11 * g14 / ( m_1 * m_1 - m2 ) + g21 * g24 / ( m_2 * m_2 - m2 ) +
543 g31 * g34 / ( m_3 * m_3 - m2 ) + g41 * g44 / ( m_4 * m_4 - m2 ) +
544 g51 * g54 / ( m_5 * m_5 - m2 ) + fs14 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
545 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
546 double_t km15 = ( g11 * g15 / ( m_1 * m_1 - m2 ) + g21 * g25 / ( m_2 * m_2 - m2 ) +
547 g31 * g35 / ( m_3 * m_3 - m2 ) + g41 * g45 / ( m_4 * m_4 - m2 ) +
548 g51 * g55 / ( m_5 * m_5 - m2 ) + fs15 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
549 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
550
551 double_t km21 = km12, km31 = km13, km41 = km14, km51 = km15;
552 double_t km23 = ( g12 * g13 / ( m_1 * m_1 - m2 ) + g22 * g23 / ( m_2 * m_2 - m2 ) +
553 g32 * g33 / ( m_3 * m_3 - m2 ) + g42 * g43 / ( m_4 * m_4 - m2 ) +
554 g52 * g53 / ( m_5 * m_5 - m2 ) ) *
555 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
556 double_t km24 = ( g12 * g14 / ( m_1 * m_1 - m2 ) + g22 * g24 / ( m_2 * m_2 - m2 ) +
557 g32 * g34 / ( m_3 * m_3 - m2 ) + g42 * g44 / ( m_4 * m_4 - m2 ) +
558 g52 * g54 / ( m_5 * m_5 - m2 ) ) *
559 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
560 double_t km25 = ( g12 * g15 / ( m_1 * m_1 - m2 ) + g22 * g25 / ( m_2 * m_2 - m2 ) +
561 g32 * g35 / ( m_3 * m_3 - m2 ) + g42 * g45 / ( m_4 * m_4 - m2 ) +
562 g52 * g55 / ( m_5 * m_5 - m2 ) ) *
563 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
564 double_t km32 = km23, km42 = km24, km52 = km25;
565 double_t km34 = ( g13 * g14 / ( m_1 * m_1 - m2 ) + g23 * g24 / ( m_2 * m_2 - m2 ) +
566 g33 * g34 / ( m_3 * m_3 - m2 ) + g43 * g44 / ( m_4 * m_4 - m2 ) +
567 g53 * g54 / ( m_5 * m_5 - m2 ) ) *
568 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
569 double_t km35 = ( g13 * g15 / ( m_1 * m_1 - m2 ) + g23 * g25 / ( m_2 * m_2 - m2 ) +
570 g33 * g35 / ( m_3 * m_3 - m2 ) + g43 * g45 / ( m_4 * m_4 - m2 ) +
571 g53 * g55 / ( m_5 * m_5 - m2 ) ) *
572 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
573 double_t km43 = km34, km53 = km35;
574 double_t km45 = ( g14 * g15 / ( m_1 * m_1 - m2 ) + g24 * g25 / ( m_2 * m_2 - m2 ) +
575 g34 * g35 / ( m_3 * m_3 - m2 ) + g44 * g45 / ( m_4 * m_4 - m2 ) +
576 g54 * g55 / ( m_5 * m_5 - m2 ) ) *
577 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
578 double_t km54 = km45;
579 double_t km22 = ( g12 * g12 / ( m_1 * m_1 - m2 ) + g22 * g22 / ( m_2 * m_2 - m2 ) +
580 g32 * g32 / ( m_3 * m_3 - m2 ) + g42 * g42 / ( m_4 * m_4 - m2 ) +
581 g52 * g52 / ( m_5 * m_5 - m2 ) ) *
582 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
583 double_t km33 = ( g13 * g13 / ( m_1 * m_1 - m2 ) + g23 * g23 / ( m_2 * m_2 - m2 ) +
584 g33 * g33 / ( m_3 * m_3 - m2 ) + g43 * g43 / ( m_4 * m_4 - m2 ) +
585 g53 * g53 / ( m_5 * m_5 - m2 ) ) *
586 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
587 double_t km44 = ( g14 * g14 / ( m_1 * m_1 - m2 ) + g24 * g24 / ( m_2 * m_2 - m2 ) +
588 g34 * g34 / ( m_3 * m_3 - m2 ) + g44 * g44 / ( m_4 * m_4 - m2 ) +
589 g54 * g54 / ( m_5 * m_5 - m2 ) ) *
590 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
591 double_t km55 = ( g15 * g15 / ( m_1 * m_1 - m2 ) + g25 * g25 / ( m_2 * m_2 - m2 ) +
592 g35 * g35 / ( m_3 * m_3 - m2 ) + g45 * g45 / ( m_4 * m_4 - m2 ) +
593 g55 * g55 / ( m_5 * m_5 - m2 ) ) *
594 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
595
596 TComplex fp11( 1.41620 * cos( 5.90540 ), 1.41620 * sin( 5.90540 ) );
597 TComplex fp12( 15.3640 * cos( 2.22510 ), 15.3640 * sin( 2.22510 ) );
598 TComplex fp13( 8.2297 * cos( 2.30550 ), 8.2297 * sin( 2.30550 ) );
599 TComplex fp14( 0.0, 0.0 );
600 TComplex fp15( 0.0, 0.0 );
601
602 TComplex beta1( 3.70140 * cos( 5.0818 ), 3.70140 * sin( 5.0818 ) );
603 TComplex beta2( 4.04470 * cos( 0.13605 ), 4.04470 * sin( 0.13605 ) );
604 TComplex beta3( 3.4591 * cos( -5.6845 ), 3.4591 * sin( -5.6845 ) );
605 TComplex beta4( 0.0, 0.0 );
606 TComplex beta5( 0.0, 0.0 );
607
608 TComplex P1 = fp11 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g11 / ( m_1 * m_1 - m2 ) +
609 beta2 * g21 / ( m_2 * m_2 - m2 ) + beta3 * g31 / ( m_3 * m_3 - m2 ) +
610 beta4 * g41 / ( m_4 * m_4 - m2 ) + beta5 * g51 / ( m_5 * m_5 - m2 );
611 TComplex P2 = fp12 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g12 / ( m_1 * m_1 - m2 ) +
612 beta2 * g22 / ( m_2 * m_2 - m2 ) + beta3 * g32 / ( m_3 * m_3 - m2 ) +
613 beta4 * g42 / ( m_4 * m_4 - m2 ) + beta5 * g52 / ( m_5 * m_5 - m2 );
614 TComplex P3 = fp13 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g13 / ( m_1 * m_1 - m2 ) +
615 beta2 * g23 / ( m_2 * m_2 - m2 ) + beta3 * g33 / ( m_3 * m_3 - m2 ) +
616 beta4 * g43 / ( m_4 * m_4 - m2 ) + beta5 * g53 / ( m_5 * m_5 - m2 );
617 TComplex P4 = fp14 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g14 / ( m_1 * m_1 - m2 ) +
618 beta2 * g24 / ( m_2 * m_2 - m2 ) + beta3 * g34 / ( m_3 * m_3 - m2 ) +
619 beta4 * g44 / ( m_4 * m_4 - m2 ) + beta5 * g54 / ( m_5 * m_5 - m2 );
620 TComplex P5 = fp15 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g15 / ( m_1 * m_1 - m2 ) +
621 beta2 * g25 / ( m_2 * m_2 - m2 ) + beta3 * g35 / ( m_3 * m_3 - m2 ) +
622 beta4 * g45 / ( m_4 * m_4 - m2 ) + beta5 * g55 / ( m_5 * m_5 - m2 );
623
624 TMatrixD MI( 5, 5 ), MA( 5, 5 ), MA_invt( 5, 5 ), MB( 5, 5 ), KM( 5, 5 ), RhoA( 5, 5 ),
625 RhoB( 5, 5 ), MRe( 5, 5 ), MIm( 5, 5 );
626 MI.UnitMatrix();
627 RhoA.UnitMatrix();
628 RhoB.UnitMatrix();
629 TMatrixDRow( KM, 0 )( 0 ) = km11;
630 TMatrixDRow( KM, 1 )( 0 ) = km21;
631 TMatrixDRow( KM, 2 )( 0 ) = km31;
632 TMatrixDRow( KM, 3 )( 0 ) = km41;
633 TMatrixDRow( KM, 4 )( 0 ) = km51;
634 TMatrixDRow( KM, 0 )( 1 ) = km12;
635 TMatrixDRow( KM, 1 )( 1 ) = km22;
636 TMatrixDRow( KM, 2 )( 1 ) = km32;
637 TMatrixDRow( KM, 3 )( 1 ) = km42;
638 TMatrixDRow( KM, 4 )( 1 ) = km52;
639 TMatrixDRow( KM, 0 )( 2 ) = km13;
640 TMatrixDRow( KM, 1 )( 2 ) = km23;
641 TMatrixDRow( KM, 2 )( 2 ) = km33;
642 TMatrixDRow( KM, 3 )( 2 ) = km43;
643 TMatrixDRow( KM, 4 )( 2 ) = km53;
644 TMatrixDRow( KM, 0 )( 3 ) = km14;
645 TMatrixDRow( KM, 1 )( 3 ) = km24;
646 TMatrixDRow( KM, 2 )( 3 ) = km34;
647 TMatrixDRow( KM, 3 )( 3 ) = km44;
648 TMatrixDRow( KM, 4 )( 3 ) = km54;
649 TMatrixDRow( KM, 0 )( 4 ) = km15;
650 TMatrixDRow( KM, 1 )( 4 ) = km25;
651 TMatrixDRow( KM, 2 )( 4 ) = km35;
652 TMatrixDRow( KM, 3 )( 4 ) = km45;
653 TMatrixDRow( KM, 4 )( 4 ) = km55;
654
655 if ( ( 1.0 - 4.0 * mass_Pion * mass_Pion / m2 ) > 0 )
656 {
657 TMatrixDRow( RhoA, 0 )( 0 ) = sqrt( 1.0 - 4.0 * mass_Pion * mass_Pion / m2 );
658 TMatrixDRow( RhoB, 0 )( 0 ) = 0.0;
659 }
660 else
661 {
662 TMatrixDRow( RhoA, 0 )( 0 ) = 0.0;
663 TMatrixDRow( RhoB, 0 )( 0 ) = sqrt( 4.0 * mass_Pion * mass_Pion / m2 - 1.0 );
664 }
665
666 if ( ( 1.0 - 4.0 * mass_Kaon * mass_Kaon / m2 ) > 0 )
667 {
668 TMatrixDRow( RhoA, 1 )( 1 ) = sqrt( 1.0 - 4.0 * mass_Kaon * mass_Kaon / m2 );
669 TMatrixDRow( RhoB, 1 )( 1 ) = 0.0;
670 }
671 else
672 {
673 TMatrixDRow( RhoA, 1 )( 1 ) = 0.0;
674 TMatrixDRow( RhoB, 1 )( 1 ) = sqrt( 4.0 * mass_Kaon * mass_Kaon / m2 - 1.0 );
675 }
676
677 TMatrixDRow( RhoA, 2 )( 2 ) = CalRho4pi( m2 );
678 TMatrixDRow( RhoB, 2 )( 2 ) = 0.0;
679 if ( ( 1.0 - 4.0 * mass_Eta * mass_Eta / m2 ) > 0 )
680 {
681 TMatrixDRow( RhoA, 3 )( 3 ) = sqrt( 1.0 - 4.0 * mass_Eta * mass_Eta / m2 );
682 TMatrixDRow( RhoB, 3 )( 3 ) = 0.0;
683 }
684 else
685 {
686 TMatrixDRow( RhoA, 3 )( 3 ) = 0.0;
687 TMatrixDRow( RhoB, 3 )( 3 ) = sqrt( 4.0 * mass_Eta * mass_Eta / m2 - 1.0 );
688 }
689
690 if ( ( 1.0 - ( mass_Eta + mass_Etap ) * ( mass_Eta + mass_Etap ) / m2 ) > 0 )
691 {
692 TMatrixDRow( RhoA, 4 )( 4 ) =
693 sqrt( 1.0 - ( mass_Eta + mass_Etap ) * ( mass_Eta + mass_Etap ) / m2 );
694 TMatrixDRow( RhoB, 4 )( 4 ) = 0.0;
695 }
696 else
697 {
698 TMatrixDRow( RhoA, 4 )( 4 ) = 0.0;
699 TMatrixDRow( RhoB, 4 )( 4 ) =
700 sqrt( ( mass_Eta + mass_Etap ) * ( mass_Eta + mass_Etap ) / m2 - 1.0 );
701 }
702
703 MA = MI + KM * RhoB;
704 MB = -1.0 * KM * RhoA;
705
706 MA_invt = MA;
707 MA_invt.Invert();
708
709 MRe = MA + MB * MA_invt * MB;
710 MRe.Invert();
711 MIm = MA_invt * MB * MRe;
712
713 TComplex ciR( 1.0, 0.0 );
714 TComplex ciM( 0.0, 1.0 );
715 TComplex amp;
716
717 amp = ( ciR * TMatrixDRow( MRe, 0 )( 0 ) - ciM * TMatrixDRow( MIm, 0 )( 0 ) ) * P1 +
718 ( ciR * TMatrixDRow( MRe, 0 )( 1 ) - ciM * TMatrixDRow( MIm, 0 )( 1 ) ) * P2 +
719 ( ciR * TMatrixDRow( MRe, 0 )( 2 ) - ciM * TMatrixDRow( MIm, 0 )( 2 ) ) * P3 +
720 ( ciR * TMatrixDRow( MRe, 0 )( 3 ) - ciM * TMatrixDRow( MIm, 0 )( 3 ) ) * P4 +
721 ( ciR * TMatrixDRow( MRe, 0 )( 4 ) - ciM * TMatrixDRow( MIm, 0 )( 4 ) ) * P5;
722
723 return amp;
724}
725
726void EvtDTopipi0pi0::calEva( double* Pic, double* Pi01, double* Pi02, double* mass1,
727 double* width1, double* amp, double* phase, int* g0, int* spin,
728 int* modetype, int nstates, double& Result, double* r0,
729 double* r1 ) {
730 double P12[4], P23[4], P13[4];
731 // double Pic[4], Pi01[4], Pi02[4], P12[4], P23[4], P13[4];
732 double cof[2], amp_PDF[2], PDF[2];
733 double spi01, spi02, scpi;
734 double Fv[2];
735 double s12, s13, s23;
736 TComplex tmpswave;
737 double Realpipis, Imagpipis;
738 for ( int i = 0; i < 4; i++ )
739 {
740 P23[i] = Pi01[i] + Pi02[i];
741 P12[i] = Pic[i] + Pi01[i];
742 P13[i] = Pic[i] + Pi02[i];
743 }
744 scpi = SCADot( Pic, Pic );
745 spi01 = SCADot( Pi01, Pi01 );
746 spi02 = SCADot( Pi02, Pi02 );
747 s23 = SCADot( P23, P23 );
748 s12 = SCADot( P12, P12 );
749 s13 = SCADot( P13, P13 ); // pipi0
750
751 double spi012[4] = { mass_Pion * mass_Pion, spi02, mass_Kaon * mass_Kaon, mk0 * mk0 };
752
753 double pro[2], temp_PDF, pro1[2], temp_PDF1, pro2[2], temp_PDF2, amp_tmp[2];
754 double mass1sq;
755 amp_PDF[0] = 0;
756 amp_PDF[1] = 0;
757 PDF[0] = 0;
758 PDF[1] = 0;
759 amp_tmp[0] = 0;
760 amp_tmp[1] = 0;
761 for ( int i = 0; i < nstates; i++ )
762 {
763 amp_tmp[0] = 0;
764 amp_tmp[1] = 0;
765 mass1sq = mass1[i] * mass1[i];
766 cof[0] = amp[i] * cos( phase[i] );
767 cof[1] = amp[i] * sin( phase[i] );
768 temp_PDF = 0;
769 if ( modetype[i] == 23 )
770 { // pi0pi0
771 temp_PDF = DDalitz( Pi01, Pi02, Pic, spin[i], mass1[i], r0[i], r1[i] );
772 if ( g0[i] == 2 ) propagatorsigma500( s23, spi01, spi02, pro );
773 if ( g0[i] == 1 )
774 propagatorRBW( mass1sq, mass1[i], width1[i], s23, spi01, spi02, r0[i] * r0[i], spin[i],
775 pro );
776 if ( g0[i] == 3 ) propagator980( mass1[i], s23, spi012, pro );
777 if ( g0[i] == 6 )
778 {
779 tmpswave = ResonanceSkm( s23 );
780 Realpipis = tmpswave.Re();
781 Imagpipis = tmpswave.Im();
782 amp_tmp[0] = temp_PDF * Realpipis;
783 amp_tmp[1] = temp_PDF * Imagpipis;
784 }
785 if ( g0[i] == 0 )
786 {
787 pro[0] = 1;
788 pro[1] = 0;
789 }
790 if ( g0[i] != 6 ) amp_tmp[0] = temp_PDF * pro[0];
791 if ( g0[i] != 6 ) amp_tmp[1] = temp_PDF * pro[1];
792 }
793 if ( modetype[i] == 12 )
794 { // pipi0
795 temp_PDF1 = DDalitz( Pic, Pi01, Pi02, spin[i], mass1[i], r0[i], r1[i] );
796 if ( g0[i] == 1 )
797 propagatorRBW( mass1sq, mass1[i], width1[i], s12, scpi, spi01, r0[i] * r0[i], spin[i],
798 pro1 );
799 if ( g0[i] == 2 )
800 propagatorGS( mass1sq, mass1[i], width1[i], s12, scpi, spi01, r0[i] * r0[i], pro1 );
801 if ( g0[i] == 0 )
802 {
803 pro1[0] = 1;
804 pro1[1] = 0;
805 }
806
807 temp_PDF2 = DDalitz( Pic, Pi02, Pi01, spin[i], mass1[i], r0[i], r1[i] );
808 if ( g0[i] == 1 )
809 propagatorRBW( mass1sq, mass1[i], width1[i], s13, scpi, spi02, r0[i] * r0[i], spin[i],
810 pro2 );
811 if ( g0[i] == 2 )
812 propagatorGS( mass1sq, mass1[i], width1[i], s13, scpi, spi02, r0[i] * r0[i], pro2 );
813 if ( g0[i] == 0 )
814 {
815 pro2[0] = 1;
816 pro2[1] = 0;
817 }
818 amp_tmp[0] = temp_PDF1 * pro1[0] + temp_PDF2 * pro2[0];
819
820 amp_tmp[1] = temp_PDF1 * pro1[1] + temp_PDF2 * pro2[1];
821 }
822 Com_Multi( amp_tmp, cof, amp_PDF );
823 PDF[0] += amp_PDF[0];
824 PDF[1] += amp_PDF[1];
825 }
826 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
827 if ( value <= 0 ) value = 1e-20;
828 Result = value;
829}
double mass
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
complex< double > TComplex
Definition Dalitz.h:14
TF1 * g1
EvtComplex exp(const EvtComplex &c)
*******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
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
virtual ~EvtDTopipi0pi0()
void decay(EvtParticle *p)
EvtDecayBase * clone()
void getName(std::string &name)
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
void checkSpinParent(EvtSpinType::spintype sp)
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)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
double get(int i) const
double double * m2
Definition qcdloop1.h:83
int t()
Definition t.c:1