BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKpipi0.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: EvtD0ToKpipi0.cc
11// the necessary file: EvtD0ToKpipi0.hh
12//
13// Description: D0 -> K- pi+ pi0, see
14// https://indico.ihep.ac.cn/event/21489/contributions/151925/attachments/76712/95007/kpipi0.pdf
15//
16// Modification history:
17//
18// Liaoyuan Dong Sat Apr 29 03:10:55 CST 2023
19// Sun May 26 17:23:25 CST 2024
20//------------------------------------------------------------------------
21#include "EvtD0ToKpipi0.hh"
22#include "../EvtGenBase/EvtComplex.hh"
23#include "../EvtGenBase/EvtDecayTable.hh"
24#include "../EvtGenBase/EvtGenKine.hh"
25#include "../EvtGenBase/EvtPDL.hh"
26#include "../EvtGenBase/EvtParticle.hh"
27#include "../EvtGenBase/EvtPatches.hh"
28#include "../EvtGenBase/EvtReport.hh"
29#include <cmath>
30#include <iostream>
31#include <stdlib.h>
32using namespace std;
33
34EvtD0ToKpipi0::~EvtD0ToKpipi0() {}
35
36void EvtD0ToKpipi0::getName( std::string& model_name ) { model_name = "D0ToKpipi0"; }
37
38EvtDecayBase* EvtD0ToKpipi0::clone() { return new EvtD0ToKpipi0; }
39
40void EvtD0ToKpipi0::init() {
41 checkNArg( 0 );
42 checkNDaug( 3 );
43 checkSpinParent( EvtSpinType::SCALAR );
44
45 mrho_770 = 0.76519, mrho_1450 = 1.465, mKn892 = 0.89555, mKp892 = 0.89167, mK0_1430 = 1.466,
46 mK2_1430 = 1.432, mKn_1680 = 1.717;
47 Grho_770 = 0.150, Grho_1450 = 0.310, GKn892 = 0.05, GKp892 = 0.0514, GK0_1430 = 0.174,
48 GK2_1430 = 0.109, GKn_1680 = 0.322;
49
50 rho[0] = 1; // rho770-pipi0
51 rho[1] = 1.5481; // rho1450-pipi0
52 rho[2] = 0.41186; // K*892p-Kpi
53 rho[3] = 0.38950; // K*982n-Kpi0
54 rho[4] = 0.37221; // K0*1430n-Kpi
55 rho[5] = 0.31480; // K0*1430p-Kpi0
56 rho[6] = 0.35958; // K2*1430n-Kpi
57 rho[7] = 0.36683; // K2*1430p-Kpi0
58 rho[8] = 1.0967; // K*1680n-Kpi
59 rho[9] = 0.26035; // K*1680p-Kpi
60 rho[10] = 1.0210; // res-S
61 rho[11] = 0.91605; // res-Kpi0p
62
63 phi[0] = 0; // rho770
64 phi[1] = 8.8326; // rho1450-pipi0
65 phi[2] = 2.9453; // K*892p-Kpi
66 phi[3] = 6.4385; // K*982n-Kpi0
67 phi[4] = 0.42526; // K0*1430n-Kpi
68 phi[5] = -1.8798; // K0*1430p-Kpi0
69 phi[6] = 5.7080; // K2*1430n-Kpi
70 phi[7] = 2.6428; // K2*1430p-Kpi0
71 phi[8] = 1.0485; // K*1680n-Kpi
72 phi[9] = 1.2300; // K*1680p-Kpi
73 phi[10] = 10.459; // res-S
74 phi[11] = 10.449; // res-Kpi0p
75
76 spin[0] = 1;
77 spin[1] = 1;
78 spin[2] = 1;
79 spin[3] = 1;
80 spin[4] = 0;
81 spin[5] = 0;
82 spin[6] = 2;
83 spin[7] = 2;
84 spin[8] = 1;
85 spin[9] = 1;
86 spin[10] = 0;
87 spin[11] = 1;
88
89 modetype[0] = 23;
90 modetype[1] = 23;
91 modetype[2] = 13;
92 modetype[3] = 12;
93 modetype[4] = 12;
94 modetype[5] = 13;
95 modetype[6] = 12;
96 modetype[7] = 13;
97 modetype[8] = 12;
98 modetype[9] = 13;
99 modetype[10] = 12;
100 modetype[11] = 13;
101
102 /*
103 std::cout << "EvtD0ToKpipi0 (May 01, 2023) ==> Initialization" << std::endl;
104 for (int i=0; i<13; i++) {
105 cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
106 }
107 */
108 mass_Ks = 0.4977;
109 mass_Eta = 0.547862;
110 mD = 1.864;
111 rD = 5;
112 metap = 0.95778;
113 mkstr = 0.89594;
114 mk0 = 0.497611;
115 mass_Kaon = 0.493677;
116 mass_Pion = 0.13957;
117 mass_Pi0 = 0.1349768;
118 math_pi = 3.1415926;
119 pi = 3.1415926;
120 mpi = 0.13957;
121
122 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
123 for ( int i = 0; i < 4; i++ )
124 {
125 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
126 }
127}
128
129void EvtD0ToKpipi0::initProbMax() {
130 setProbMax( 188.35 ); // setProbMax;
131}
132
133void EvtD0ToKpipi0::decay( EvtParticle* p ) {
134 //-----------for max value------------------
135 /*
136 double maxprob = 0.0;
137 for(int ir=0;ir<=60000000;ir++){
138 p->initializePhaseSpace(getNDaug(),getDaugs());
139 EvtVector4R ks0 = p->getDaug(0)->getP4();
140 EvtVector4R pic = p->getDaug(1)->getP4();
141 EvtVector4R pi0 = p->getDaug(2)->getP4();
142
143 double Ks0[4],Pic[4],Pi0[4];
144 Ks0[0] = ks0.get(0); Pic[0] = pic.get(0); Pi0[0] = pi0.get(0);
145 Ks0[1] = ks0.get(1); Pic[1] = pic.get(1); Pi0[1] = pi0.get(1);
146 Ks0[2] = ks0.get(2); Pic[2] = pic.get(2); Pi0[2] = pi0.get(2);
147 Ks0[3] = ks0.get(3); Pic[3] = pic.get(3); Pi0[3] = pi0.get(3);
148 double value;
149 calPDF(Ks0, Pic, Pi0, value);
150 if(value>maxprob) {
151 maxprob=value;
152 std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
153 }
154 }
155 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
156 return;
157 */
158 //-----------------------------------------------
159 p->initializePhaseSpace( getNDaug(), getDaugs() );
160 EvtVector4R ks0 = p->getDaug( 0 )->getP4();
161 EvtVector4R pic = p->getDaug( 1 )->getP4();
162 EvtVector4R pi0 = p->getDaug( 2 )->getP4();
163
164 double Ks0[4], Pic[4], Pi0[4];
165 Ks0[0] = ks0.get( 0 );
166 Pic[0] = pic.get( 0 );
167 Pi0[0] = pi0.get( 0 );
168 Ks0[1] = ks0.get( 1 );
169 Pic[1] = pic.get( 1 );
170 Pi0[1] = pi0.get( 1 );
171 Ks0[2] = ks0.get( 2 );
172 Pic[2] = pic.get( 2 );
173 Pi0[2] = pi0.get( 2 );
174 Ks0[3] = ks0.get( 3 );
175 Pic[3] = pic.get( 3 );
176 Pi0[3] = pi0.get( 3 );
177
178 double value;
179 calPDF( Ks0, Pic, Pi0, value );
180
181 setProb( value );
182 return;
183}
184
185double EvtD0ToKpipi0::calPDF( double Ks0[], double Pic[], double Pi0[], double& Result ) {
186 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
187 double flag[3], mass_R[2], width_R[2];
188 double pro[2];
189 double sa[3], sb[3], sc[3], B[3];
190 double t1D[4], t1V[4], t1A[4];
191 double pS[4], pV[4], pD[4], pA[4];
192 double pro1[2], pro2[2], proKPi_S[2];
193
194 double rD2 = 25.0;
195 double rRes2 = 9.0;
196 double rRes = 9.0;
197 double mass1[12] = { mrho_770, mrho_1450, mKp892, mKn892, mK0_1430, mK0_1430,
198 mK2_1430, mK2_1430, mKn_1680, mKn_1680, 1.414, 1.414 };
199 double width1[12] = { Grho_770, Grho_1450, GKp892, GKn892, GK0_1430, GK0_1430,
200 GK2_1430, GK2_1430, GKn_1680, GKn_1680, 0.232, 0.232 };
201 double g0[12] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0 };
202 double temp_PDF = 0;
203 double P12[4], P23[4], P13[4];
204 double scpi, snpi, sks0;
205 double s12, s13, s23;
206 for ( int i = 0; i < 4; i++ )
207 {
208 P12[i] = Pic[i] + Ks0[i];
209 P13[i] = Pi0[i] + Ks0[i];
210 P23[i] = Pic[i] + Pi0[i];
211 }
212 scpi = SCADot( Pic, Pic );
213 snpi = SCADot( Pi0, Pi0 );
214 sks0 = SCADot( Ks0, Ks0 );
215 s12 = SCADot( P12, P12 );
216 s13 = SCADot( P13, P13 );
217 s23 = SCADot( P23, P23 );
218 double mass1sq;
219 double Amp_KPiS[2];
220 amp_PDF[0] = 0;
221 amp_PDF[1] = 0;
222 PDF[0] = 0;
223 PDF[1] = 0;
224
225 for ( int i = 0; i < 12; i++ )
226 {
227 amp_tmp[0] = 0;
228 amp_tmp[1] = 0;
229 mass1sq = mass1[i] * mass1[i];
230 cof[0] = rho[i] * cos( phi[i] );
231 cof[1] = rho[i] * sin( phi[i] );
232 temp_PDF = 0;
233 if ( modetype[i] == 23 )
234 {
235 temp_PDF = DDalitz( Pic, Pi0, Ks0, spin[i], mass1[i] );
236 if ( g0[i] == 1 )
237 propagatorGS( mass1sq, mass1[i], width1[i], s23, scpi, snpi, 9.0, pro );
238 // if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],s23,scpi,snpi,rRes,spin[i],pro);
239 if ( g0[i] == 2 ) KPiSLASS( s23, scpi, snpi, pro );
240 if ( g0[i] == 3 ) propagatorFlatte( mass1[i], width1[i], s23, scpi, snpi, 0, pro );
241 if ( g0[i] == 0 )
242 {
243 pro[0] = 1;
244 pro[1] = 0;
245 }
246 amp_tmp[0] = temp_PDF * pro[0];
247 amp_tmp[1] = temp_PDF * pro[1];
248 }
249 if ( modetype[i] == 12 )
250 {
251 temp_PDF = DDalitz( Ks0, Pic, Pi0, spin[i], mass1[i] );
252 if ( g0[i] == 1 )
253 propagatorRBW( mass1sq, mass1[i], width1[i], s12, sks0, scpi, rRes, spin[i], pro );
254 if ( g0[i] == 2 )
255 {
256 propagatorFlatte( mass1[i], width1[i], s12, sks0, scpi, 1, pro );
257 pro[0] = pro[0] * ( 0.01 + 0.990 * 0.990 ) / ( 0.01 + s12 );
258 pro[1] = pro[1] * ( 0.01 + 0.990 * 0.990 ) / ( 0.01 + s12 );
259 }
260 if ( g0[i] == 3 )
261 propagatorFlatte( mass1[i], width1[i], s12, sks0, scpi, 1, pro ); // Only for a0(980)
262
263 if ( g0[i] == 0 )
264 {
265 pro[0] = 1;
266 pro[1] = 0;
267 }
268 amp_tmp[0] = temp_PDF * pro[0];
269 amp_tmp[1] = temp_PDF * pro[1];
270 }
271 if ( modetype[i] == 13 )
272 {
273 temp_PDF = DDalitz( Ks0, Pi0, Pic, spin[i], mass1[i] );
274 if ( g0[i] == 1 )
275 propagatorRBW( mass1sq, mass1[i], width1[i], s13, sks0, snpi, rRes, spin[i], pro );
276 if ( g0[i] == 2 )
277 { // K*1430 Flatte
278 double skm2[2] = { sks0, mass_Ks * mass_Ks };
279 double spi2[2] = { snpi, mass_Eta * mass_Eta };
280 propagatorKstr1430( mass1[i], s13, skm2, spi2, pro );
281 };
282 if ( g0[i] == 3 )
283 propagatorFlatte( mass1[i], width1[i], s13, sks0, snpi, 1, pro ); // Only for a0(980)
284 if ( g0[i] == 0 )
285 {
286 pro[0] = 1;
287 pro[1] = 0;
288 }
289 amp_tmp[0] = temp_PDF * pro[0];
290 amp_tmp[1] = temp_PDF * pro[1];
291 }
292 if ( modetype[i] == 132 )
293 {
294 KPiSLASS( s12, sks0, scpi, Amp_KPiS );
295 amp_tmp[0] = Amp_KPiS[0];
296 amp_tmp[1] = Amp_KPiS[1];
297 }
298 Com_Multi( amp_tmp, cof, amp_PDF );
299 PDF[0] += amp_PDF[0];
300 PDF[1] += amp_PDF[1];
301 }
302 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
303 if ( value <= 0 ) value = 1e-20;
304 Result = value;
305
306 return value;
307}
308
309void EvtD0ToKpipi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
310 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
311 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
312}
313void EvtD0ToKpipi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
314 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
315 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
316 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
317}
318//------------base---------------------------------
319double EvtD0ToKpipi0::SCADot( double a1[4], double a2[4] ) {
320 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
321 return _cal;
322}
323
324double EvtD0ToKpipi0::barrier( int l, double sa, double sb, double sc, double r,
325 double mass ) {
326 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
327 if ( q < 0 ) q = 1e-16;
328 double z;
329 z = q * r * r;
330 double sa0;
331 sa0 = mass * mass;
332 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
333 if ( q0 < 0 ) q0 = 1e-16;
334 double z0 = q0 * r * r;
335 double F = 0.0;
336 if ( l == 0 ) F = 1;
337 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
338 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
339 return F;
340}
341
342double EvtD0ToKpipi0::Barrier( int l, double sa, double sb, double sc, double r2 ) {
343 double F;
344 double tmp = sa + sb - sc;
345 double q = 0.25 * tmp * tmp / sa - sb;
346 if ( q < 0 ) q = 1e-16;
347 double z = q * r2;
348 if ( l == 1 ) { F = sqrt( 2.0 * z / ( 1.0 + z ) ); }
349 else if ( l == 2 )
350 {
351 double z2 = z * z;
352 F = sqrt( 13.0 * z2 / ( 9.0 + 3.0 * z + z2 ) );
353 }
354 else { F = 1.0; }
355 return F;
356}
357
358//------------------spin-------------------------------------------
359void EvtD0ToKpipi0::calt1( double daug1[4], double daug2[4], double t1[4] ) {
360 double p, pq, tmp;
361 double pa[4], qa[4];
362 for ( int i = 0; i < 4; i++ )
363 {
364 pa[i] = daug1[i] + daug2[i];
365 qa[i] = daug1[i] - daug2[i];
366 }
367 p = SCADot( pa, pa );
368 pq = SCADot( pa, qa );
369 tmp = pq / p;
370 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
371}
372void EvtD0ToKpipi0::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
373 double p, r;
374 double pa[4], t1[4];
375 calt1( daug1, daug2, t1 );
376 r = SCADot( t1, t1 ) / 3.0;
377 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
378 p = SCADot( pa, pa );
379 for ( int i = 0; i < 4; i++ )
380 {
381 for ( int j = 0; j < 4; j++ )
382 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
383 }
384}
385//-------------------prop--------------------------------------------
386void EvtD0ToKpipi0::propagator( double mass2, double mass, double width, double sx,
387 double prop[2] ) {
388 double a[2], b[2];
389 a[0] = 1;
390 a[1] = 0;
391 b[0] = mass2 - sx;
392 b[1] = -mass * width;
393 Com_Divide( a, b, prop );
394}
395double EvtD0ToKpipi0::wid( double mass2, double mass, double sa, double sb, double sc,
396 double r2, int l ) {
397 double widm = 0.;
398 double m = sqrt( sa );
399 double tmp = sb - sc;
400 double tmp1 = sa + tmp;
401 double q = 0.25 * tmp1 * tmp1 / sa - sb;
402 if ( q < 0 ) q = 1e-16;
403 double tmp2 = mass2 + tmp;
404 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
405 if ( q0 < 0 ) q0 = 1e-16;
406 double z = q * r2;
407 double z0 = q0 * r2;
408 double t = q / q0;
409 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
410 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
411 else if ( l == 2 )
412 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
413 return widm;
414}
415double EvtD0ToKpipi0::widl1( double mass2, double mass, double sa, double sb, double sc,
416 double r2 ) {
417 double widm = 0.;
418 double m = sqrt( sa );
419 double tmp = sb - sc;
420 double tmp1 = sa + tmp;
421 double q = 0.25 * tmp1 * tmp1 / sa - sb;
422 if ( q < 0 ) q = 1e-16;
423 double tmp2 = mass2 + tmp;
424 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
425 if ( q0 < 0 ) q0 = 1e-16;
426 double z = q * r2;
427 double z0 = q0 * r2;
428 double F = ( 1 + z0 ) / ( 1 + z );
429 double t = q / q0;
430 widm = t * sqrt( t ) * mass / m * F;
431 return widm;
432}
433void EvtD0ToKpipi0::propagatorRBW( double mass2, double mass, double width, double sa,
434 double sb, double sc, double r2, int l, double prop[2] ) {
435 double a[2], b[2];
436 a[0] = 1;
437 a[1] = 0;
438 double q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
439
440 b[0] = mass2 - sa;
441 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
442 Com_Divide( a, b, prop );
443}
444
445void EvtD0ToKpipi0::propagatorFlatte( double mass, double width, double sa, double sb,
446 double sc, int r, double prop[2] ) {
447 double q, qKsK, qetapi;
448 // double qKsK,qetapi;
449 double rhoab[2], rhoKsK[2];
450 q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
451 qetapi = 0.25 * ( sa + 0.547862 * 0.547862 - 0.13957039 * 0.13957039 ) *
452 ( sa + 0.547862 * 0.547862 - 0.13957039 * 0.13957039 ) / sa -
453 0.547862 * 0.547862;
454 if ( r == 0 ) qKsK = 0.25 * sa - 0.49368 * 0.49368;
455 if ( r == 1 ) qKsK = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
456 if ( qetapi > 0 )
457 {
458 rhoab[0] = 2 * sqrt( qetapi / sa );
459 rhoab[1] = 0;
460 }
461 if ( qetapi < 0 )
462 {
463 rhoab[0] = 0;
464 rhoab[1] = 2 * sqrt( -qetapi / sa );
465 }
466 if ( qKsK > 0 )
467 {
468 rhoKsK[0] = 2 * sqrt( qKsK / sa );
469 rhoKsK[1] = 0;
470 }
471 if ( qKsK < 0 )
472 {
473 rhoKsK[0] = 0;
474 rhoKsK[1] = 2 * sqrt( -qKsK / sa );
475 }
476 double a[2], b[2];
477 a[0] = 1;
478 a[1] = 0;
479 b[0] = mass * mass - sa + 0.341 * rhoab[1] + 0.892 * 0.341 * rhoKsK[1];
480 b[1] = -( 0.341 * rhoab[0] + 0.892 * 0.341 * rhoKsK[0] );
481 Com_Divide( a, b, prop );
482}
483void EvtD0ToKpipi0::PiPiSWASS( double sa, double sb, double sc, double prop[2] ) {
484 double tmp = sb - sc;
485 double tmp2 = sa + tmp;
486 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
487 double q = sqrt( qs );
488 double a0 = -0.11 / mass_Pion;
489 prop[0] = 1 / ( 1 + a0 * a0 * q * q );
490 prop[1] = a0 * q / ( 1 + a0 * a0 * q * q );
491}
492
493void EvtD0ToKpipi0::propagatorGS( double mass2, double mass, double width, double sa,
494 double sb, double sc, double r2, double prop[2] ) {
495 double GS1 = 0.636619783;
496 double GS2 = 0.01860182466;
497 double GS3 = 0.1591549458; // 1/(2*math_2pi)
498 double GS4 = 0.00620060822; // mass_Pion2/math_pi
499 double a[2], b[2];
500 double tmp = sb - sc;
501 double tmp1 = sa + tmp;
502 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
503 if ( q2 < 0 ) q2 = 1e-16;
504
505 double tmp2 = mass2 + tmp;
506 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
507 if ( q02 < 0 ) q02 = 1e-16;
508
509 double q = sqrt( q2 );
510 double q0 = sqrt( q02 );
511 double m = sqrt( sa );
512 double q03 = q0 * q02;
513 // double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
514 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
515
516 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
517 double h0 = GS1 * q0 / mass * tmp3;
518 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
519 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
520 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
521
522 a[0] = 1.0 + d * width / mass;
523 a[1] = 0.0;
524 b[0] = mass2 - sa + width * f;
525 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
526 Com_Divide( a, b, prop );
527}
528
529void EvtD0ToKpipi0::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
530 const double m1430 = 1.441;
531 const double sa0 = 1.441 * 1.441; // m1430*m1430;
532 const double w1430 = 0.193;
533 const double Lass1 = 0.25 / sa0;
534 double tmp = sb - sc;
535 double tmp1 = sa0 + tmp;
536 double q0 = Lass1 * tmp1 * tmp1 - sb;
537 if ( q0 < 0 ) q0 = 1e-16;
538 double tmp2 = sa + tmp;
539 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
540 double q = sqrt( qs );
541 double width = w1430 * q * m1430 / sqrt( sa * q0 );
542 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
543 if ( temp_R < 0 ) temp_R += math_pi;
544 double deltaR = -1.915 + temp_R; // fiR=-109.7
545 double temp_F =
546 atan( 0.226 * q / ( 2.0 - 3.819 * qs ) ); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
547 if ( temp_F < 0 ) temp_F += math_pi;
548 double deltaF = 0.002 + temp_F; // fiF=0.1
549 double deltaS = deltaR + 2.0 * deltaF;
550 double t1 = 0.96 * sin( deltaF );
551 double t2 = sin( deltaR );
552 double CF[2], CS[2];
553 CF[0] = cos( deltaF );
554 CF[1] = sin( deltaF );
555 CS[0] = cos( deltaS );
556 CS[1] = sin( deltaS );
557 prop[0] = t1 * CF[0] + t2 * CS[0];
558 prop[1] = t1 * CF[1] + t2 * CS[1];
559}
560
561void EvtD0ToKpipi0::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
562 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
563 if ( q > 0 )
564 {
565 rho[0] = 2 * sqrt( q / sa );
566 rho[1] = 0;
567 }
568 else if ( q < 0 )
569 {
570 rho[0] = 0;
571 rho[1] = 2 * sqrt( -q / sa );
572 }
573}
574
575void EvtD0ToKpipi0::propagatorKstr1430( double mass, double sx, double* sb, double* sc,
576 double prop[2] ) // K*1430 Flatte
577{
578 double unit[2] = { 1.0 };
579 double ci[2] = { 0, 1 };
580 double rho1[2];
581 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
582 double rho2[2];
583 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
584 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
585 double tmp1[2] = { gKPi_Kstr1430, 0 };
586 double tmp11[2];
587 double tmp2[2] = { gEtaPK_Kstr1430, 0 };
588 double tmp22[2];
589 Com_Multi( tmp1, rho1, tmp11 );
590 Com_Multi( tmp2, rho2, tmp22 );
591 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
592 double tmp31[2];
593 Com_Multi( tmp3, ci, tmp31 );
594 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
595 Com_Divide( unit, tmp4, prop );
596}
597double EvtD0ToKpipi0::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
598 double mass ) {
599 double pR[4], pD[4];
600 double temp_PDF, v_re;
601 temp_PDF = 0.0;
602 v_re = 0.0;
603 double B[2], s1, s2, s3, sR, sD;
604 for ( int i = 0; i < 4; i++ )
605 {
606 pR[i] = P1[i] + P2[i];
607 pD[i] = pR[i] + P3[i];
608 }
609 s1 = SCADot( P1, P1 );
610 s2 = SCADot( P2, P2 );
611 s3 = SCADot( P3, P3 );
612 sR = SCADot( pR, pR );
613 sD = SCADot( pD, pD );
614 int G[4][4];
615 for ( int i = 0; i != 4; i++ )
616 {
617 for ( int j = 0; j != 4; j++ )
618 {
619 if ( i == j )
620 {
621 if ( i == 0 ) G[i][j] = 1;
622 else G[i][j] = -1;
623 }
624 else G[i][j] = 0;
625 }
626 }
627 if ( Ang == 0 )
628 {
629 B[0] = 1;
630 B[1] = 1;
631 temp_PDF = 1;
632 }
633 if ( Ang == 1 )
634 {
635 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
636 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.864 );
637 double t1[4], T1[4];
638 calt1( P1, P2, t1 );
639 calt1( pR, P3, T1 );
640 temp_PDF = 0;
641 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
642 }
643 if ( Ang == 2 )
644 {
645 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
646 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.864 );
647 double t2[4][4], T2[4][4];
648 calt2( P1, P2, t2 );
649 calt2( pR, P3, T2 );
650 temp_PDF = 0;
651 for ( int i = 0; i < 4; i++ )
652 {
653 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
654 }
655 }
656 v_re = temp_PDF * B[0] * B[1];
657 return v_re;
658}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
a0
character *LEPTONflag integer iresonances real zeta5 real a0
Definition Ekhara/src/ekhara/common.h:4
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
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
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
EvtD0ToKpipi0::~EvtD0ToKpipi0
virtual ~EvtD0ToKpipi0()
Definition EvtD0ToKpipi0.cc:34
EvtD0ToKpipi0::clone
EvtDecayBase * clone()
Definition EvtD0ToKpipi0.cc:38
EvtD0ToKpipi0::decay
void decay(EvtParticle *p)
Definition EvtD0ToKpipi0.cc:133
EvtD0ToKpipi0::getName
void getName(std::string &name)
Definition EvtD0ToKpipi0.cc:36
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:479
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:261
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:465
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:68
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:449
EvtDecayProb::setProb
void setProb(double prob)
Definition EvtDecayProb.hh:40
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:108
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1042
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:171
t
int t()
Definition t.c:1