BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKKpipi0.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: EvtDToKKpipi0.cc
11// the necessary file: EvtDToKKpipi0.hh
12//
13// Description: D+ -> K- K+ pi+ pi0
14//
15// Modification history:
16// Liaoyuan Dong, Zhenxuan Li May 26 16:55:36 2024 Module created
17//------------------------------------------------------------------------
18#include "EvtDToKKpipi0.hh"
19#include "../EvtGenBase/EvtComplex.hh"
20#include "../EvtGenBase/EvtDecayTable.hh"
21#include "../EvtGenBase/EvtGenKine.hh"
22#include "../EvtGenBase/EvtPDL.hh"
23#include "../EvtGenBase/EvtParticle.hh"
24#include "../EvtGenBase/EvtPatches.hh"
25#include "../EvtGenBase/EvtReport.hh"
26#include "TComplex.h"
27#include "TMath.h"
28#include "TMatrixD.h"
29#include "TSpline.h"
30#include <cmath>
31#include <iostream>
32#include <stdio.h>
33#include <stdlib.h>
34using namespace std;
35
36EvtDToKKpipi0::~EvtDToKKpipi0() {}
37
38void EvtDToKKpipi0::getName( std::string& model_name ) { model_name = "DToKKpipi0"; }
39
40EvtDecayBase* EvtDToKKpipi0::clone() { return new EvtDToKKpipi0; }
41
42void EvtDToKKpipi0::init() {
43
44 checkNArg( 0 );
45 checkNDaug( 4 );
46 checkSpinParent( EvtSpinType::SCALAR );
47
48 mD = 1.86966; // D+ mass
49 // mD = 1.86483; //D0 mass
50 mDs = 1.9683;
51 rRes = 3.0;
52 rD = 5.0;
53 metap = 0.95778;
54 mkstr = 0.89594;
55 mk0 = 0.497614;
56 mass_Kaon = 0.49368;
57 mass_Pion = 0.13957;
58 mass_Eta = 0.547862;
59 mass_Etap = 0.95778;
60 mass_Pi0 = 0.1349766;
61 math_pi = 3.1415926;
62 mrho = 0.77511;
63 Grho = 1.2037e-01;
64 mrho1450 = 1.465;
65 Grho1450 = 0.400;
66 ma0 = 0.990;
67 Ga0 = 0.0756;
68 ma01450 = 1.474;
69 Ga01450 = 0.265;
70 meta = 0.547862;
71 ma2 = 1.3177;
72 ma20 = 1.698;
73 Ga2 = 0.1111;
74 Ga20 = 0.265;
75 mf980 = 0.965;
76 Gf980 = 0.0756;
77 meta1405 = 1.4088;
78 Geta1405 = 0.0501;
79 //---------------------------------
80 mphi = 1.019461;
81 Gphi = 0.004249;
82 // Gphi = 6.6912e-03;
83 mnKstr = 0.89555;
84 GnKstr = 0.0473;
85 mcKstr = 0.89167;
86 GcKstr = 0.0514;
87 mK1270 = 1.28981;
88 mK1400 = 1.403;
89 mf1285 = 1.2819;
90 Gf1285 = 0.0227;
91 mf1420 = 1.4263;
92 Gf1420 = 0.0545;
93 gpieta2 = 0.341;
94 gKK2 = 0.892 * 0.341;
95
96 phi[0] = 0.0;
97 rho[0] = 1.0; // D2VV 12340 phi rho S
98 phi[1] = -2.4183895412;
99 rho[1] = 0.3660444327; // D2VV 12340 phi rho P
100 phi[2] = 3.9126761898;
101 rho[2] = 0.9905707322; // D2VV 12340 phi rho D
102 phi[3] = -3.4581051622;
103 rho[3] = 0.4970069076;
104 phi[4] = 0.6059414627;
105 rho[4] = 0.4744002304;
106 phi[5] = 0.7501566868;
107 rho[5] = 0.6407798868;
108 phi[6] = -1.3470518401;
109 rho[6] = -2.2701523251;
110 phi[7] = -0.9948321319;
111 rho[7] = -1.0841075539;
112 phi[8] = -0.7144319026;
113 rho[8] = 1.3650576749;
114 phi[9] = -0.0298815830;
115 rho[9] = -0.3101468271;
116 phi[10] = -1.0962771498;
117 rho[10] = 0.3799755719;
118 phi[11] = 0.2459805626;
119 rho[11] = -0.3419295418;
120 phi[12] = -1.4937738501;
121 rho[12] = -0.4275950103;
122
123 sp1270 = new TSpline3( "sp1270", s, width1270, 1206 );
124 sp1400 = new TSpline3( "sp1400", s, width1400, 1206 );
125
126 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
127 int EE[4][4][4][4] = {
128 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
129 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
130 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
131 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
132 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
133 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
134 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
135 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
136 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
137 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
138 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
139 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
140 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
141 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
142 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
143 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
144 for ( int i = 0; i < 4; i++ )
145 {
146 for ( int j = 0; j < 4; j++ )
147 {
148 G[i][j] = GG[i][j];
149 for ( int k = 0; k < 4; k++ )
150 {
151 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
152 }
153 }
154 }
155}
156
157void EvtDToKKpipi0::initProbMax() { setProbMax( 301000 ); }
158
159void EvtDToKKpipi0::decay( EvtParticle* p ) {
160 //-----------for max value------------------
161 /*double maxprob = 0.0;
162 for(int ir=0;ir<=60000000;ir++){
163 p->initializePhaseSpace(getNDaug(),getDaugs());
164 EvtVector4R p1 = p->getDaug(0)->getP4();
165 EvtVector4R p2 = p->getDaug(1)->getP4();
166 EvtVector4R p3 = p->getDaug(2)->getP4();
167 EvtVector4R p4 = p->getDaug(3)->getP4();
168 double value;
169 double P1[4],P2[4],P3[4],P4[4];
170 mother_c=EvtPDL::getStdHep(p->getId());
171 if(mother_c==411){
172 P1[0] = p1.get(0); P2[0] = p2.get(0); P3[0] = p3.get(0); P4[0] = p4.get(0);
173 P1[1] = p1.get(1); P2[1] = p2.get(1); P3[1] = p3.get(1); P4[1] = p4.get(1);
174 P1[2] = p1.get(2); P2[2] = p2.get(2); P3[2] = p3.get(2); P4[2] = p4.get(2);
175 P1[3] = p1.get(3); P2[3] = p2.get(3); P3[3] = p3.get(3); P4[3] = p4.get(3);
176 }else if(mother_c==-411){
177 P1[0] = p1.get(0); P2[0] = p2.get(0); P3[0] = p3.get(0); P4[0] = p4.get(0);
178 P1[1] = -1.0*p1.get(1); P2[1] = -1.0*p2.get(1); P3[1] = -1.0*p3.get(1); P4[1]
179 = -1.0*p4.get(1); P1[2] = -1.0*p1.get(2); P2[2] = -1.0*p2.get(2); P3[2] = -1.0*p3.get(2);
180 P4[2] = -1.0*p4.get(2); P1[3] = -1.0*p1.get(3); P2[3] = -1.0*p2.get(3); P3[3] =
181 -1.0*p3.get(3); P4[3] = -1.0*p4.get(3);
182 }
183 calPDF(P1, P2, P3, P4, value);
184 if(value>maxprob) {
185 maxprob=value;
186 std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
187 }
188 }
189 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
190 return;*/
191 //-----------------------------------------------
192 p->initializePhaseSpace( getNDaug(), getDaugs() );
193 EvtVector4R p1 = p->getDaug( 0 )->getP4();
194 EvtVector4R p2 = p->getDaug( 1 )->getP4();
195 EvtVector4R p3 = p->getDaug( 2 )->getP4();
196 EvtVector4R p4 = p->getDaug( 3 )->getP4();
197 mother_c = EvtPDL::getStdHep( p->getId() );
198
199 double P1[4], P2[4], P3[4], P4[4];
200 if ( mother_c == 411 )
201 {
202 P1[0] = p1.get( 0 );
203 P2[0] = p2.get( 0 );
204 P3[0] = p3.get( 0 );
205 P4[0] = p4.get( 0 );
206 P1[1] = p1.get( 1 );
207 P2[1] = p2.get( 1 );
208 P3[1] = p3.get( 1 );
209 P4[1] = p4.get( 1 );
210 P1[2] = p1.get( 2 );
211 P2[2] = p2.get( 2 );
212 P3[2] = p3.get( 2 );
213 P4[2] = p4.get( 2 );
214 P1[3] = p1.get( 3 );
215 P2[3] = p2.get( 3 );
216 P3[3] = p3.get( 3 );
217 P4[3] = p4.get( 3 );
218 }
219 else if ( mother_c == -411 )
220 {
221 P1[0] = p1.get( 0 );
222 P2[0] = p2.get( 0 );
223 P3[0] = p3.get( 0 );
224 P4[0] = p4.get( 0 );
225 P1[1] = -1.0 * p1.get( 1 );
226 P2[1] = -1.0 * p2.get( 1 );
227 P3[1] = -1.0 * p3.get( 1 );
228 P4[1] = -1.0 * p4.get( 1 );
229 P1[2] = -1.0 * p1.get( 2 );
230 P2[2] = -1.0 * p2.get( 2 );
231 P3[2] = -1.0 * p3.get( 2 );
232 P4[2] = -1.0 * p4.get( 2 );
233 P1[3] = -1.0 * p1.get( 3 );
234 P2[3] = -1.0 * p2.get( 3 );
235 P3[3] = -1.0 * p3.get( 3 );
236 P4[3] = -1.0 * p4.get( 3 );
237 }
238
239 double value;
240 calPDF( P1, P2, P3, P4, value );
241 setProb( value );
242 return;
243}
244
245void EvtDToKKpipi0::calPDF( double P1[], double P2[], double P3[], double P4[],
246 double& Result ) {
247
248 double s12, s13, s14, s23, s24, s34;
249 double s123, s124, s134, s234;
250
251 double P12[4], P13[4], P14[4], P23[4], P24[4], P34[4];
252 double P123[4], P124[4], P134[4], P234[4];
253 for ( int i = 0; i < 4; i++ )
254 {
255 P12[i] = P1[i] + P2[i];
256 P13[i] = P1[i] + P3[i];
257 P14[i] = P1[i] + P4[i];
258 P23[i] = P2[i] + P3[i];
259 P24[i] = P2[i] + P4[i];
260 P34[i] = P3[i] + P4[i];
261 P123[i] = P12[i] + P3[i];
262 P124[i] = P12[i] + P4[i];
263 P134[i] = P13[i] + P4[i];
264 P234[i] = P23[i] + P4[i];
265 }
266
267 s12 = dot( P12, P12 );
268 s13 = dot( P13, P13 );
269 s14 = dot( P14, P14 );
270 s23 = dot( P23, P23 );
271 s24 = dot( P24, P24 );
272 s34 = dot( P34, P34 );
273
274 s123 = dot( P123, P123 );
275 s124 = dot( P124, P124 );
276 s134 = dot( P134, P134 );
277 s234 = dot( P234, P234 );
278 TComplex PDF[100], cof, pdf, module, pro1, pro2, cpro1, cpro2;
279 //---------
280 double wid1270, wid1400;
281 wid1270 = sp1270->Eval( s134 );
282 wid1400 = sp1400->Eval( s134 );
283 //---------
284
285 pro1 = propagatorRBW( mphi, Gphi, s12, sck, sck, 3.0, 1, 1 );
286 pro2 = propagatorGS( mrho, Grho, s34, scpi, snpi, 3.0, 1 );
287 PDF[0] = Spin_factor( P1, P2, P3, P4, 0, mphi, mrho, 3.0, 5.0, 0 ) * pro1 *
288 pro2; // phi rho S-wave
289 PDF[1] = Spin_factor( P1, P2, P3, P4, 1, mphi, mrho, 3.0, 5.0, 0 ) * pro1 *
290 pro2; // phi rho P-wave
291 PDF[2] = Spin_factor( P1, P2, P3, P4, 2, mphi, mrho, 3.0, 5.0, 0 ) * pro1 *
292 pro2; // phi rho D-wave
293
294 pro1 = propagatorRBW( mnKstr, GnKstr, s13, sck, scpi, 3.0, 1, 1 );
295 pro2 = propagatorRBW( mcKstr, GcKstr, s24, sck, snpi, 3.0, 1, 1 );
296 PDF[3] = Spin_factor( P1, P3, P2, P4, 0, mnKstr, mcKstr, 3.0, 5.0, 0 ) * pro1 *
297 pro2; // K*0bar K*+ S-wave
298 PDF[4] = Spin_factor( P1, P3, P2, P4, 1, mnKstr, mcKstr, 3.0, 5.0, 0 ) * pro1 *
299 pro2; // K*0bar K*+ P-wave
300 PDF[5] = Spin_factor( P1, P3, P2, P4, 2, mnKstr, mcKstr, 3.0, 5.0, 0 ) * pro1 *
301 pro2; // K*0bar K*+ D-wave
302
303 pro1 = propagatorGS( mrho, Grho, s34, scpi, snpi, 3.0, 1 );
304 pro2 = propagatorRBW( mK1270, wid1270, s134, s34, sck, 3.0, 1, 0 );
305 PDF[6] = Spin_factor( P3, P4, P1, P2, 0, mrho, mK1270, 3.0, 5.0, 1 ) * pro1 *
306 pro2; // K1270 pi, K1270->rho K S-wave
307
308 pro1 = propagatorRBW( mnKstr, GnKstr, s13, sck, scpi, 3.0, 1, 1 );
309 pro2 = propagatorRBW( mK1400, wid1400, s134, s13, snpi, 3.0, 1, 0 );
310 cpro1 = propagatorRBW( mcKstr, GcKstr, s14, sck, snpi, 3.0, 1, 1 );
311 cpro2 = propagatorRBW( mK1400, wid1400, s134, s14, scpi, 3.0, 1, 0 );
312 PDF[7] = Spin_factor( P1, P3, P4, P2, 0, mnKstr, mK1400, 3.0, 5.0, 1 ) * pro1 * pro2 +
313 Spin_factor( P1, P4, P3, P2, 0, mcKstr, mK1400, 3.0, 5.0, 1 ) * cpro1 *
314 cpro2; // K1400 pi, K1400->Kstr pi
315
316 pro1 = propagatorGS( mrho, Grho, s34, scpi, snpi, 3.0, 1 );
317 double sb[2], sc[2], pro_tmp[2];
318 // sb[0] = snpi; sc[0] = seta;
319 // sb[1] = sck; sc[1] = sck;
320 // pro2 = propagatorFlatte(ma0, Ga0, s12, sb, sc, gpieta2, gKK2);
321 propagatorFlatte( ma0, Ga0, s12, snpi, seta, 0, pro_tmp );
322 pro2 = TComplex( pro_tmp[0], pro_tmp[1] );
323 PDF[8] = Spin_factor( P3, P4, P1, P2, 0, mrho, ma0, 3.0, 5.0, 10 ) * pro1 * pro2; // rho a0
324
325 pro1 = propagatorRBW( mcKstr, GcKstr, s14, sck, snpi, 3.0, 1, 1 );
326 pro2 = propagatorRBW( mf1420, Gf1420, s124, s14, sck, 3.0, 0, 1 );
327 cpro1 = propagatorRBW( mcKstr, GcKstr, s24, sck, snpi, 3.0, 1, 1 );
328 cpro2 = propagatorRBW( mf1420, Gf1420, s124, s24, sck, 3.0, 0, 1 );
329 PDF[9] = Spin_factor( P1, P4, P2, P3, 0, mcKstr, mf1420, 3.0, 5.0, 1 ) * pro1 * pro2 +
330 Spin_factor( P2, P4, P1, P3, 0, mcKstr, mf1420, 3.0, 5.0, 1 ) * cpro1 *
331 cpro2; // f1420 -> Kstr K
332
333 pro1 = propagatorRBW( mcKstr, GcKstr, s14, sck, snpi, 3.0, 1, 1 );
334 pro2 = propagatorRBW( mf1285, Gf1285, s124, s14, sck, 3.0, 0, 1 );
335 cpro1 = propagatorRBW( mcKstr, GcKstr, s24, sck, snpi, 3.0, 1, 1 );
336 cpro2 = propagatorRBW( mf1285, Gf1285, s124, s24, sck, 3.0, 0, 1 );
337 PDF[10] = Spin_factor( P1, P4, P2, P3, 0, mcKstr, mf1285, 3.0, 5.0, 1 ) * pro1 * pro2 +
338 Spin_factor( P2, P4, P1, P3, 0, mcKstr, mf1285, 3.0, 5.0, 1 ) * cpro1 *
339 cpro2; // f1285 -> Kstr K
340
341 // pro1 = propagatorFlatte(ma0, Ga0, s12, sb, sc, gpieta2, gKK2);
342 propagatorFlatte( ma0, Ga0, s12, snpi, seta, 0, pro_tmp );
343 pro1 = TComplex( pro_tmp[0], pro_tmp[1] );
344 pro2 = propagatorRBW( meta1405, Geta1405, s124, s12, snpi, 3.0, 1, 1 );
345 PDF[11] = pro1 * pro2; // eta1405->a0 pi
346
347 pro1 = propagatorRBW( mcKstr, GcKstr, s14, sck, snpi, 3.0, 1, 1 );
348 pro2 = TComplex( pro_tmp[0], pro_tmp[1] );
349 pro2 = propagatorRBW( meta1405, Geta1405, s124, s14, sck, 3.0, 1, 1 );
350 cpro1 = propagatorRBW( mcKstr, GcKstr, s24, sck, snpi, 3.0, 1, 1 );
351 cpro2 = propagatorRBW( meta1405, Geta1405, s124, s24, sck, 3.0, 1, 1 );
352 PDF[12] = Spin_factor( P1, P4, P2, P3, 1, mcKstr, meta1405, 3.0, 5.0, 21 ) * pro1 * pro2 +
353 Spin_factor( P2, P4, P1, P3, 1, mcKstr, meta1405, 3.0, 5.0, 21 ) * cpro1 *
354 cpro2; // eta1405 -> Kstr K
355
356 pdf = TComplex( 0, 0 );
357 for ( Int_t i = 0; i < 13; i++ )
358 {
359 cof = TComplex( rho[i] * TMath::Cos( phi[i] ), rho[i] * TMath::Sin( phi[i] ) );
360 // cout<<i<<" "<<(cof*PDF[i]).Re()<<" "<<(cof*PDF[i]).Im()<<endl;
361 pdf += cof * PDF[i];
362 }
363 module = pdf * TComplex::Conjugate( pdf );
364 double value;
365 value = module.Re();
366 // cout<<value<<endl;
367 // return (value <= 0) ? 1e-20 : value;
368 Result = ( value <= 0 ) ? 1e-20 : value;
369}
370
371TComplex EvtDToKKpipi0::propogator( double mass, double width, double sx ) const {
372 TComplex ci( 0, 1 );
373 TComplex prop = 1.0 / ( mass * mass - sx - ci * mass * width );
374 return prop;
375}
376double EvtDToKKpipi0::wid( double mass, double sa, double sb, double sc, double r,
377 int l ) const {
378 double widm( 0. ), q( 0. ), q0( 0. );
379 double sa0 = mass * mass;
380 double m = sqrt( sa );
381 q = Qabcs( sa, sb, sc );
382 q0 = Qabcs( sa0, sb, sc );
383 double z, z0;
384 z = q * r * r;
385 z0 = q0 * r * r;
386 double t = q / q0;
387 double F( 0. );
388 if ( l == 0 ) F = 1;
389 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
390 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
391 widm = pow( t, l + 0.5 ) * mass / m * F * F;
392 return widm;
393}
394double EvtDToKKpipi0::h( double m, double q ) const {
395 double h( 0. );
396 h = 2 / math_pi * q / m * log( ( m + 2 * q ) / ( 0.13957 + 0.134976 ) );
397 return h;
398}
399double EvtDToKKpipi0::dh( double mass, double q0 ) const {
400 double dh = h( mass, q0 ) * ( 1.0 / ( 8 * q0 * q0 ) - 1.0 / ( 2 * mass * mass ) ) +
401 1.0 / ( 2 * math_pi * mass * mass );
402 return dh;
403}
404double EvtDToKKpipi0::f( double mass, double sx, double q0, double q ) const {
405 double m = sqrt( sx );
406 double f = mass * mass / ( pow( q0, 3 ) ) *
407 ( q * q * ( h( m, q ) - h( mass, q0 ) ) +
408 ( mass * mass - sx ) * q0 * q0 * dh( mass, q0 ) );
409 return f;
410}
411double EvtDToKKpipi0::d( double mass, double q0 ) const {
412 double cmpi = 0.5 * ( 0.13957 + 0.134976 );
413 // double cmpi = 0.13957;
414 double d =
415 3.0 / math_pi * cmpi * cmpi / ( q0 * q0 ) * log( ( mass + 2 * q0 ) / ( 2 * cmpi ) ) +
416 mass / ( 2 * math_pi * q0 ) - ( cmpi * cmpi * mass ) / ( math_pi * pow( q0, 3 ) );
417 return d;
418}
419TComplex EvtDToKKpipi0::propagatorRBW( double mass, double width, double sa, double sb,
420 double sc, double r, int l, int isvaried ) const {
421 TComplex ci( 0, 1 );
422 TComplex prop;
423 if ( isvaried == 1 )
424 prop = 1.0 / ( mass * mass - sa - ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
425 if ( isvaried == 0 ) prop = 1.0 / ( mass * mass - sa - ci * mass * width );
426 return prop;
427}
428TComplex EvtDToKKpipi0::propagatorGS( double mass, double width, double sa, double sb,
429 double sc, double r, int l ) const {
430 TComplex ci( 0, 1 );
431 double q = Qabcs( sa, sb, sc );
432 double sa0 = mass * mass;
433 double q0 = Qabcs( sa0, sb, sc );
434 q = sqrt( q );
435 q0 = sqrt( q0 );
436 TComplex prop = ( 1 + d( mass, q0 ) * width / mass ) /
437 ( mass * mass - sa + width * f( mass, sa, q0, q ) -
438 ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
439 // TComplex prop =
440 // (1+d(mass,q0)*width/mass)/(mass*mass-sa+width*f(mass,sa,q0,q)-ci*mass*width);
441 return prop;
442}
443TComplex EvtDToKKpipi0::Flatte_rhoab( double sa, double sb, double sc ) const {
444 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
445 TComplex rho;
446 TComplex ci( 0, 1 );
447 if ( q > 0 ) { rho = 2.0 * sqrt( q / sa ) * ( TComplex::One() ); }
448 if ( q < 0 ) { rho = 2.0 * sqrt( -q / sa ) * ci; }
449 return rho;
450}
451void EvtDToKKpipi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
452 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
453 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
454}
455void EvtDToKKpipi0::propagatorFlatte( double mass, double width, double sa, double sb,
456 double sc, int r, double prop[2] ) {
457 // r == 0, neutral a0
458 // r == 1, charged a0
459 double q, qKsK;
460 double rhoab[2], rhoKsK[2];
461 q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
462 if ( r == 0 ) qKsK = 0.25 * sa - 0.49368 * 0.49368;
463 if ( r == 1 )
464 qKsK =
465 0.25 * ( sa + 3.899750596e-03 ) * ( sa + 3.899750596e-03 ) / sa - 0.497614 * 0.497614;
466 if ( q > 0 )
467 {
468 rhoab[0] = 2 * sqrt( q / sa );
469 rhoab[1] = 0;
470 }
471 if ( q < 0 )
472 {
473 rhoab[0] = 0;
474 rhoab[1] = 2 * sqrt( -q / sa );
475 }
476 if ( qKsK > 0 )
477 {
478 rhoKsK[0] = 2 * sqrt( qKsK / sa );
479 rhoKsK[1] = 0;
480 }
481 if ( qKsK < 0 )
482 {
483 rhoKsK[0] = 0;
484 rhoKsK[1] = 2 * sqrt( -qKsK / sa );
485 }
486 double a[2], b[2];
487 a[0] = 1;
488 a[1] = 0;
489 b[0] = mass * mass - sa + 0.341 * rhoab[1] + 0.892 * 0.341 * rhoKsK[1];
490 b[1] = -( 0.341 * rhoab[0] + 0.892 * 0.341 * rhoKsK[0] );
491 Com_Divide( a, b, prop );
492}
493double EvtDToKKpipi0::Spin_factor( double P1[], double P2[], double P3[], double P4[],
494 Int_t Ang, double mass1, double mass2, double rRes,
495 double rD, Int_t mode ) {
496 // mode == 0, D->VV
497 // mode == 1, D->AP, A->VP, V->PP
498 // mode == 2, D->VP, V->VP, V->PP
499
500 // mode == 10, D->VS, V->PP, S->PP
501 // mode == 11, D->TS, T->PP, S->PP
502 // mode == 12, PHSP
503 // mode == 20, D->AP, A->SP, S->PP
504 // mode == 21, D->P'P, P'->VP, V->PP
505
506 double pR1[4], pR2[4], pD[4];
507 double temp_PDF, v_re;
508 temp_PDF = 0.0;
509 v_re = 0.0;
510 double mD = 1.86966;
511 double B[3], s1, s2, s3, s4, sR1, sR2, sD;
512 sD = mD * mD;
513 double t1D[4], t2D[4][4];
514 s1 = dot( P1, P1 );
515 s2 = dot( P2, P2 );
516 s3 = dot( P3, P3 );
517 s4 = dot( P4, P4 );
518 double tV1[4], tV2[4];
519
520 //=============================
521 if ( mode == 0 )
522 {
523 // VV
524 for ( int i = 0; i < 4; i++ )
525 {
526 pR1[i] = P1[i] + P2[i];
527 pR2[i] = P3[i] + P4[i];
528 pD[i] = pR1[i] + pR2[i];
529 }
530 sR1 = dot( pR1, pR1 );
531 sR2 = dot( pR2, pR2 );
532 calt1( P1, P2, tV1 );
533 calt1( P3, P4, tV2 );
534 calt1( pR1, pR2, t1D );
535 calt2( pR1, pR2, t2D );
536 B[0] = barrier( 1, sR1, s1, s2, rRes, mass1 );
537 B[1] = barrier( 1, sR2, s3, s4, rRes, mass2 );
538
539 if ( Ang == 0 )
540 {
541 B[2] = 1.0;
542 // for(int i=0; i<4; i++){
543 // temp_PDF += tV1[i]*tV2[i]*G[i][i];
544 // }
545 for ( int i = 0; i < 4; i++ )
546 {
547 for ( int j = 0; j < 4; j++ )
548 {
549 temp_PDF += ( -G[i][j] + pD[i] * pD[j] / sD ) * tV1[i] * tV2[j] * G[i][i] * G[j][j];
550 }
551 }
552 }
553 if ( Ang == 1 )
554 {
555 B[2] = barrier( 1, sD, sR1, sR2, rD, mD );
556 for ( int i = 0; i < 4; i++ )
557 {
558 for ( int j = 0; j < 4; j++ )
559 {
560 for ( int k = 0; k < 4; k++ )
561 {
562 for ( int l = 0; l < 4; l++ )
563 {
564 temp_PDF += E[i][j][k][l] * pD[i] * t1D[j] * tV1[k] * tV2[l] * G[i][i] *
565 G[j][j] * G[k][k] * G[l][l];
566 }
567 }
568 }
569 }
570 }
571 if ( Ang == 2 )
572 {
573 B[2] = barrier( 2, sD, sR1, sR2, rD, mD );
574 for ( int i = 0; i < 4; i++ )
575 {
576 for ( int j = 0; j < 4; j++ )
577 {
578 // temp_PDF +=
579 // pD[i]*pD[j]*tV1[i]*tV2[j]*G[i][i]*G[j][j];
580 temp_PDF += t2D[i][j] * tV1[i] * tV2[j] * G[i][i] * G[j][j];
581 }
582 }
583 }
584 }
585 //=============================
586 if ( mode == 1 )
587 {
588 // AP
589 for ( int i = 0; i < 4; i++ )
590 {
591 pR1[i] = P1[i] + P2[i];
592 pR2[i] = pR1[i] + P3[i];
593 pD[i] = pR2[i] + P4[i];
594 }
595 double tV[4];
596 sR1 = dot( pR1, pR1 );
597 sR2 = dot( pR2, pR2 );
598 calt1( P1, P2, tV );
599 calt1( pR2, P4, t1D );
600 B[0] = barrier( 1, sR1, s1, s2, rRes, mass1 );
601 B[2] = barrier( 1, sD, sR2, s4, rD, mD );
602 if ( Ang == 0 )
603 {
604 B[1] = 1.0;
605 for ( int i = 0; i < 4; i++ )
606 {
607 for ( int j = 0; j < 4; j++ )
608 {
609 temp_PDF +=
610 t1D[i] * ( -G[i][j] + pR2[i] * pR2[j] / sR2 ) * tV[j] * G[i][i] * G[j][j];
611 }
612 }
613 }
614 if ( Ang == 2 )
615 {
616 B[1] = barrier( 2, sR2, sR1, s3, rRes, mass2 );
617 double t2A[4][4];
618 calt2( pR1, P3, t2A );
619 for ( int i = 0; i < 4; i++ )
620 {
621 for ( int j = 0; j < 4; j++ )
622 { temp_PDF += t1D[i] * t2A[i][j] * tV[j] * G[i][i] * G[j][j]; }
623 }
624 }
625 }
626 //=============================
627 if ( mode == 2 )
628 {
629 // VP
630 for ( int i = 0; i < 4; i++ )
631 {
632 pR1[i] = P1[i] + P2[i];
633 pR2[i] = pR1[i] + P3[i];
634 pD[i] = pR2[i] + P4[i];
635 }
636 sR1 = dot( pR1, pR1 );
637 sR2 = dot( pR2, pR2 );
638 B[0] = barrier( 1, sR1, s1, s2, rRes, mass1 );
639 B[1] = barrier( 1, sR2, sR1, s3, rRes, mass2 );
640 B[2] = barrier( 1, sD, sR2, s4, rD, mD );
641 for ( int i = 0; i < 4; i++ )
642 {
643 for ( int j = 0; j < 4; j++ )
644 {
645 for ( int k = 0; k < 4; k++ )
646 {
647 for ( int l = 0; l < 4; l++ )
648 {
649 temp_PDF += E[i][j][k][l] * pR2[i] * ( pR1[j] - P3[j] ) * P4[k] *
650 ( P1[l] - P2[l] ) * G[i][i] * G[j][j] * G[k][k] * G[l][l];
651 }
652 }
653 }
654 }
655 }
656 //=============================
657 if ( mode == 10 )
658 {
659 // VS
660 for ( int i = 0; i < 4; i++ )
661 {
662 pR1[i] = P1[i] + P2[i];
663 pR2[i] = P3[i] + P4[i];
664 pD[i] = pR1[i] + pR2[i];
665 }
666 sR1 = dot( pR1, pR1 );
667 sR2 = dot( pR2, pR2 );
668 calt1( P1, P2, tV1 );
669 calt1( pR1, pR2, t1D );
670 B[0] = barrier( 1, sR1, s1, s2, rRes, mass1 );
671 B[1] = 1.0;
672 B[2] = barrier( 1, sD, sR1, sR2, rD, mD );
673 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1D[i] * tV1[i] * G[i][i]; }
674 }
675 //=============================
676 if ( mode == 11 )
677 {
678 // TS
679 for ( int i = 0; i < 4; i++ )
680 {
681 pR1[i] = P1[i] + P2[i];
682 pR2[i] = P3[i] + P4[i];
683 pD[i] = pR1[i] + pR2[i];
684 }
685 sR1 = dot( pR1, pR1 );
686 sR2 = dot( pR2, pR2 );
687 double tT[4][4], t2D[4][4];
688 calt2( P1, P2, tT );
689 calt2( pR1, pR2, t2D );
690 B[0] = barrier( 2, sR1, s1, s2, rRes, mass1 );
691 B[1] = 1.0;
692 B[2] = barrier( 2, sD, sR1, sR2, rD, mD );
693 for ( int i = 0; i < 4; i++ )
694 {
695 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2D[i][j] * tT[i][j] * G[i][i] * G[j][j]; }
696 }
697 }
698 //=============================
699 if ( mode == 20 )
700 {
701 // AP, A->SP
702 for ( int i = 0; i < 4; i++ )
703 {
704 pR1[i] = P1[i] + P2[i];
705 pR2[i] = pR1[i] + P3[i];
706 pD[i] = pR2[i] + P4[i];
707 }
708 double tA[4];
709 sR1 = dot( pR1, pR1 );
710 sR2 = dot( pR2, pR2 );
711 calt1( pR1, P3, tA );
712 calt1( pR2, P4, t1D );
713 B[0] = 1.0;
714 B[1] = barrier( 1, sR2, sR1, s3, rRes, mass2 );
715 B[2] = barrier( 1, sD, sR2, s4, rD, mD );
716 for ( int i = 0; i < 4; i++ ) { temp_PDF += tA[i] * t1D[i] * G[i][i]; }
717 }
718 //=============================
719 if ( mode == 21 )
720 {
721 // P' P, P'->VP, V->PP
722 for ( int i = 0; i < 4; i++ )
723 {
724 pR1[i] = P1[i] + P2[i];
725 pR2[i] = pR1[i] + P3[i];
726 pD[i] = pR2[i] + P4[i];
727 }
728 double tVP[4];
729 sR1 = dot( pR1, pR1 );
730 sR2 = dot( pR2, pR2 );
731 calt1( P1, P2, tVP );
732 B[0] = barrier( 1, sR1, s1, s2, rRes, mass1 );
733 B[1] = barrier( 1, sR2, sR1, s3, rRes, mass2 );
734 B[2] = 1.0;
735 for ( int i = 0; i < 4; i++ ) { temp_PDF += P3[i] * tVP[i] * G[i][i]; }
736 }
737 //=============================
738 v_re = temp_PDF * B[0] * B[1] * B[2];
739 return v_re;
740}
741double EvtDToKKpipi0::CalRho4pi( double s ) {
742
743 if ( s >= 1. ) { return sqrt( ( s - 16. * mass_Pion * mass_Pion ) / s ); }
744 else
745 {
746 double gam = 0.0005;
747 double s2 = s * s;
748 double s3 = s2 * s;
749 double s4 = s2 * s2;
750 double s5 = s4 * s;
751 double s6 = s5 * s;
752
753 gam -= 0.0193 * s;
754 gam += 0.1385 * s2;
755 gam -= 0.2084 * s3;
756 gam -= 0.2974 * s4;
757 gam += 0.1366 * s5;
758 gam += 1.0789 * s6;
759
760 return gam;
761 }
762}
763
764TComplex EvtDToKKpipi0::ResonanceSkm( double m2 ) {
765
766 // all par value
767 double g11 = 0.22889, g12 = -0.55377, g13 = 0.00, g14 = -0.39899, g15 = -0.34639;
768 double g21 = 0.94128, g22 = 0.55095, g23 = 0.00, g24 = 0.39065, g25 = 0.31503;
769 double g31 = 0.36856, g32 = 0.23888, g33 = 0.55639, g34 = 0.18340, g35 = 0.18681;
770 double g41 = 0.33650, g42 = 0.40907, g43 = 0.85679, g44 = 0.19906, g45 = -0.00984;
771 double g51 = 0.18171, g52 = -0.17558, g53 = -0.79658, g54 = -0.00355, g55 = 0.22358;
772
773 double fs11 = 0.23399, fs12 = 0.15044, fs13 = -0.20545, fs14 = 0.32825, fs15 = 0.35412;
774 double m_1 = 0.65100, m_2 = 1.20360, m_3 = 1.55817, m_4 = 1.21000, m_5 = 1.82206;
775 double ss0 = -3.92637;
776 double sp0 = -3.0; // v0
777 // double sp0 = -0.07;//v1
778 double sA0 = -0.15;
779 double sA = 1.0;
780
781 double km11 = ( g11 * g11 / ( m_1 * m_1 - m2 ) + g21 * g21 / ( m_2 * m_2 - m2 ) +
782 g31 * g31 / ( m_3 * m_3 - m2 ) + g41 * g41 / ( m_4 * m_4 - m2 ) +
783 g51 * g51 / ( m_5 * m_5 - m2 ) + fs11 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
784 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
785 double km12 = ( g11 * g12 / ( m_1 * m_1 - m2 ) + g21 * g22 / ( m_2 * m_2 - m2 ) +
786 g31 * g32 / ( m_3 * m_3 - m2 ) + g41 * g42 / ( m_4 * m_4 - m2 ) +
787 g51 * g52 / ( m_5 * m_5 - m2 ) + fs12 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
788 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
789 double km13 = ( g11 * g13 / ( m_1 * m_1 - m2 ) + g21 * g23 / ( m_2 * m_2 - m2 ) +
790 g31 * g33 / ( m_3 * m_3 - m2 ) + g41 * g43 / ( m_4 * m_4 - m2 ) +
791 g51 * g53 / ( m_5 * m_5 - m2 ) + fs13 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
792 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
793 double km14 = ( g11 * g14 / ( m_1 * m_1 - m2 ) + g21 * g24 / ( m_2 * m_2 - m2 ) +
794 g31 * g34 / ( m_3 * m_3 - m2 ) + g41 * g44 / ( m_4 * m_4 - m2 ) +
795 g51 * g54 / ( m_5 * m_5 - m2 ) + fs14 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
796 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
797 double km15 = ( g11 * g15 / ( m_1 * m_1 - m2 ) + g21 * g25 / ( m_2 * m_2 - m2 ) +
798 g31 * g35 / ( m_3 * m_3 - m2 ) + g41 * g45 / ( m_4 * m_4 - m2 ) +
799 g51 * g55 / ( m_5 * m_5 - m2 ) + fs15 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
800 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
801
802 double km21 = km12, km31 = km13, km41 = km14, km51 = km15;
803
804 double km23 = ( g12 * g13 / ( m_1 * m_1 - m2 ) + g22 * g23 / ( m_2 * m_2 - m2 ) +
805 g32 * g33 / ( m_3 * m_3 - m2 ) + g42 * g43 / ( m_4 * m_4 - m2 ) +
806 g52 * g53 / ( m_5 * m_5 - m2 ) ) *
807 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
808 double km24 = ( g12 * g14 / ( m_1 * m_1 - m2 ) + g22 * g24 / ( m_2 * m_2 - m2 ) +
809 g32 * g34 / ( m_3 * m_3 - m2 ) + g42 * g44 / ( m_4 * m_4 - m2 ) +
810 g52 * g54 / ( m_5 * m_5 - m2 ) ) *
811 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
812 double km25 = ( g12 * g15 / ( m_1 * m_1 - m2 ) + g22 * g25 / ( m_2 * m_2 - m2 ) +
813 g32 * g35 / ( m_3 * m_3 - m2 ) + g42 * g45 / ( m_4 * m_4 - m2 ) +
814 g52 * g55 / ( m_5 * m_5 - m2 ) ) *
815 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
816
817 double km32 = km23, km42 = km24, km52 = km25;
818 double km34 = ( g13 * g14 / ( m_1 * m_1 - m2 ) + g23 * g24 / ( m_2 * m_2 - m2 ) +
819 g33 * g34 / ( m_3 * m_3 - m2 ) + g43 * g44 / ( m_4 * m_4 - m2 ) +
820 g53 * g54 / ( m_5 * m_5 - m2 ) ) *
821 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
822 double km35 = ( g13 * g15 / ( m_1 * m_1 - m2 ) + g23 * g25 / ( m_2 * m_2 - m2 ) +
823 g33 * g35 / ( m_3 * m_3 - m2 ) + g43 * g45 / ( m_4 * m_4 - m2 ) +
824 g53 * g55 / ( m_5 * m_5 - m2 ) ) *
825 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
826 double km43 = km34, km53 = km35;
827 double km45 = ( g14 * g15 / ( m_1 * m_1 - m2 ) + g24 * g25 / ( m_2 * m_2 - m2 ) +
828 g34 * g35 / ( m_3 * m_3 - m2 ) + g44 * g45 / ( m_4 * m_4 - m2 ) +
829 g54 * g55 / ( m_5 * m_5 - m2 ) ) *
830 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
831 double km54 = km45;
832
833 double km22 = ( g12 * g12 / ( m_1 * m_1 - m2 ) + g22 * g22 / ( m_2 * m_2 - m2 ) +
834 g32 * g32 / ( m_3 * m_3 - m2 ) + g42 * g42 / ( m_4 * m_4 - m2 ) +
835 g52 * g52 / ( m_5 * m_5 - m2 ) ) *
836 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
837 double km33 = ( g13 * g13 / ( m_1 * m_1 - m2 ) + g23 * g23 / ( m_2 * m_2 - m2 ) +
838 g33 * g33 / ( m_3 * m_3 - m2 ) + g43 * g43 / ( m_4 * m_4 - m2 ) +
839 g53 * g53 / ( m_5 * m_5 - m2 ) ) *
840 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
841 double km44 = ( g14 * g14 / ( m_1 * m_1 - m2 ) + g24 * g24 / ( m_2 * m_2 - m2 ) +
842 g34 * g34 / ( m_3 * m_3 - m2 ) + g44 * g44 / ( m_4 * m_4 - m2 ) +
843 g54 * g54 / ( m_5 * m_5 - m2 ) ) *
844 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
845 double km55 = ( g15 * g15 / ( m_1 * m_1 - m2 ) + g25 * g25 / ( m_2 * m_2 - m2 ) +
846 g35 * g35 / ( m_3 * m_3 - m2 ) + g45 * g45 / ( m_4 * m_4 - m2 ) +
847 g55 * g55 / ( m_5 * m_5 - m2 ) ) *
848 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
849 TComplex fp11( 16.55, -3.87 );
850 TComplex beta1( 8.89, 3.36 );
851 TComplex fp12( -13.99, 5.12 );
852 TComplex beta2( 16.20, 5.77 );
853 TComplex fp13( -77.34, -80.04 );
854 TComplex beta3( -21.60, 27.41 );
855 TComplex fp14( -28.52, -3.20 );
856 TComplex beta4( -40.15, 35.36 );
857 TComplex fp15( 30.47, 5.70 );
858 TComplex beta5( 30.0, -43.09 );
859
860 TComplex P1 = fp11 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g11 / ( m_1 * m_1 - m2 ) +
861 beta2 * g21 / ( m_2 * m_2 - m2 ) + beta3 * g31 / ( m_3 * m_3 - m2 ) +
862 beta4 * g41 / ( m_4 * m_4 - m2 ) + beta5 * g51 / ( m_5 * m_5 - m2 );
863 TComplex P2 = fp12 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g12 / ( m_1 * m_1 - m2 ) +
864 beta2 * g22 / ( m_2 * m_2 - m2 ) + beta3 * g32 / ( m_3 * m_3 - m2 ) +
865 beta4 * g42 / ( m_4 * m_4 - m2 ) + beta5 * g52 / ( m_5 * m_5 - m2 );
866 TComplex P3 = fp13 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g13 / ( m_1 * m_1 - m2 ) +
867 beta2 * g23 / ( m_2 * m_2 - m2 ) + beta3 * g33 / ( m_3 * m_3 - m2 ) +
868 beta4 * g43 / ( m_4 * m_4 - m2 ) + beta5 * g53 / ( m_5 * m_5 - m2 );
869 TComplex P4 = fp14 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g14 / ( m_1 * m_1 - m2 ) +
870 beta2 * g24 / ( m_2 * m_2 - m2 ) + beta3 * g34 / ( m_3 * m_3 - m2 ) +
871 beta4 * g44 / ( m_4 * m_4 - m2 ) + beta5 * g54 / ( m_5 * m_5 - m2 );
872 TComplex P5 = fp15 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g15 / ( m_1 * m_1 - m2 ) +
873 beta2 * g25 / ( m_2 * m_2 - m2 ) + beta3 * g35 / ( m_3 * m_3 - m2 ) +
874 beta4 * g45 / ( m_4 * m_4 - m2 ) + beta5 * g55 / ( m_5 * m_5 - m2 );
875
876 TMatrixD MI( 5, 5 ), MA( 5, 5 ), MA_invt( 5, 5 ), MB( 5, 5 ), KM( 5, 5 ), RhoA( 5, 5 ),
877 RhoB( 5, 5 ), MRe( 5, 5 ), MIm( 5, 5 );
878 MI.UnitMatrix();
879 RhoA.UnitMatrix();
880 RhoB.UnitMatrix();
881
882 TMatrixDRow( KM, 0 )( 0 ) = km11;
883 TMatrixDRow( KM, 1 )( 0 ) = km21;
884 TMatrixDRow( KM, 2 )( 0 ) = km31;
885 TMatrixDRow( KM, 3 )( 0 ) = km41;
886 TMatrixDRow( KM, 4 )( 0 ) = km51;
887 TMatrixDRow( KM, 0 )( 1 ) = km12;
888 TMatrixDRow( KM, 1 )( 1 ) = km22;
889 TMatrixDRow( KM, 2 )( 1 ) = km32;
890 TMatrixDRow( KM, 3 )( 1 ) = km42;
891 TMatrixDRow( KM, 4 )( 1 ) = km52;
892 TMatrixDRow( KM, 0 )( 2 ) = km13;
893 TMatrixDRow( KM, 1 )( 2 ) = km23;
894 TMatrixDRow( KM, 2 )( 2 ) = km33;
895 TMatrixDRow( KM, 3 )( 2 ) = km43;
896 TMatrixDRow( KM, 4 )( 2 ) = km53;
897 TMatrixDRow( KM, 0 )( 3 ) = km14;
898 TMatrixDRow( KM, 1 )( 3 ) = km24;
899 TMatrixDRow( KM, 2 )( 3 ) = km34;
900 TMatrixDRow( KM, 3 )( 3 ) = km44;
901 TMatrixDRow( KM, 4 )( 3 ) = km54;
902 TMatrixDRow( KM, 0 )( 4 ) = km15;
903 TMatrixDRow( KM, 1 )( 4 ) = km25;
904 TMatrixDRow( KM, 2 )( 4 ) = km35;
905 TMatrixDRow( KM, 3 )( 4 ) = km45;
906 TMatrixDRow( KM, 4 )( 4 ) = km55;
907
908 TMatrixDRow( RhoA, 0 )( 0 ) = sqrt( 1.0 - 4.0 * mass_Pion * mass_Pion / m2 );
909 TMatrixDRow( RhoB, 0 )( 0 ) = 0.0;
910 if ( ( 1.0 - 4.0 * mass_Kaon * mass_Kaon / m2 ) > 0 )
911 {
912 TMatrixDRow( RhoA, 1 )( 1 ) = sqrt( 1.0 - 4.0 * mass_Kaon * mass_Kaon / m2 );
913 TMatrixDRow( RhoB, 1 )( 1 ) = 0.0;
914 }
915 else
916 {
917 TMatrixDRow( RhoA, 1 )( 1 ) = 0.0;
918 TMatrixDRow( RhoB, 1 )( 1 ) = sqrt( 4.0 * mass_Kaon * mass_Kaon / m2 - 1.0 );
919 }
920
921 TMatrixDRow( RhoA, 2 )( 2 ) = CalRho4pi( m2 );
922 TMatrixDRow( RhoB, 2 )( 2 ) = 0.0;
923 if ( ( 1.0 - 4.0 * mass_Eta * mass_Eta / m2 ) > 0 )
924 {
925 TMatrixDRow( RhoA, 3 )( 3 ) = sqrt( 1.0 - 4.0 * mass_Eta * mass_Eta / m2 );
926 TMatrixDRow( RhoB, 3 )( 3 ) = 0.0;
927 }
928 else
929 {
930 TMatrixDRow( RhoA, 3 )( 3 ) = 0.0;
931 TMatrixDRow( RhoB, 3 )( 3 ) = sqrt( 4.0 * mass_Eta * mass_Eta / m2 - 1.0 );
932 }
933 TMatrixDRow( RhoA, 4 )( 4 ) = 0.0;
934 TMatrixDRow( RhoB, 4 )( 4 ) =
935 sqrt( ( mass_Eta + mass_Etap ) * ( mass_Eta + mass_Etap ) / m2 - 1.0 );
936
937 MA = MI + KM * RhoB;
938 MB = -1.0 * KM * RhoA;
939
940 MA_invt = MA;
941 MA_invt.Invert();
942
943 MRe = MA + MB * MA_invt * MB;
944 MRe.Invert();
945 MIm = MA_invt * MB * MRe;
946
947 TComplex ciR( 1.0, 0.0 );
948 TComplex ciM( 0.0, 1.0 );
949 TComplex amp;
950
951 amp = ( ciR * TMatrixDRow( MRe, 0 )( 0 ) - ciM * TMatrixDRow( MIm, 0 )( 0 ) ) * P1 +
952 ( ciR * TMatrixDRow( MRe, 0 )( 1 ) - ciM * TMatrixDRow( MIm, 0 )( 1 ) ) * P2 +
953 ( ciR * TMatrixDRow( MRe, 0 )( 2 ) - ciM * TMatrixDRow( MIm, 0 )( 2 ) ) * P3 +
954 ( ciR * TMatrixDRow( MRe, 0 )( 3 ) - ciM * TMatrixDRow( MIm, 0 )( 3 ) ) * P4 +
955 ( ciR * TMatrixDRow( MRe, 0 )( 4 ) - ciM * TMatrixDRow( MIm, 0 )( 4 ) ) * P5;
956
957 return amp;
958}
959
960double EvtDToKKpipi0::dot( double* a1, double* a2 ) const {
961 double dot = 0;
962 for ( Int_t i = 0; i != 4; i++ ) { dot += a1[i] * a2[i] * G[i][i]; }
963 return dot;
964}
965double EvtDToKKpipi0::Qabcs( double sa, double sb, double sc ) const {
966 Double_t Qabcs = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
967 if ( Qabcs < 0 ) Qabcs = 1e-16;
968 return Qabcs;
969}
970double EvtDToKKpipi0::barrier( double l, double sa, double sb, double sc, double r,
971 double mass ) const {
972 double sa0 = mass * mass;
973 double q0 = Qabcs( sa0, sb, sc );
974 double z0 = q0 * r * r;
975 double q = Qabcs( sa, sb, sc );
976 q = sqrt( q );
977 double z = q * r;
978 z = z * z;
979 double F = 1;
980 if ( l > 2 ) F = 0;
981 if ( l == 0 ) F = 1;
982 if ( l == 1 ) { F = sqrt( ( 1 + z0 ) / ( 1 + z ) ); }
983 if ( l == 2 ) { F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) ); }
984 return F;
985}
986void EvtDToKKpipi0::calt1( double daug1[], double daug2[], double t1[] ) const {
987 double p, pq;
988 double pa[4], qa[4];
989 for ( int i = 0; i != 4; i++ )
990 {
991 pa[i] = daug1[i] + daug2[i];
992 qa[i] = daug1[i] - daug2[i];
993 }
994 p = dot( pa, pa );
995 pq = dot( pa, qa );
996 for ( int i = 0; i != 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
997}
998void EvtDToKKpipi0::calt2( Double_t daug1[], Double_t daug2[], Double_t t2[][4] ) const {
999 double p, r;
1000 double pa[4], t1[4];
1001 calt1( daug1, daug2, t1 );
1002 r = dot( t1, t1 );
1003 for ( int i = 0; i != 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
1004 p = dot( pa, pa );
1005 for ( int i = 0; i != 4; i++ )
1006 {
1007 for ( int j = 0; j != 4; j++ )
1008 { t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * ( G[i][j] - pa[i] * pa[j] / p ); }
1009 }
1010}
p2
double p2[4]
Definition BabayagaNLO.cxx:41
p1
double p1[4]
Definition BabayagaNLO.cxx:41
mass
double mass
Definition CosmicGenerator.cxx:128
TComplex
complex< double > TComplex
Definition Dalitz.h:14
s
XmlRpcServer s
Definition HelloServer.cpp:11
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
EvtDToKKpipi0::~EvtDToKKpipi0
virtual ~EvtDToKKpipi0()
Definition EvtDToKKpipi0.cc:36
EvtDToKKpipi0::clone
EvtDecayBase * clone()
Definition EvtDToKKpipi0.cc:40
EvtDToKKpipi0::getName
void getName(std::string &name)
Definition EvtDToKKpipi0.cc:38
EvtDToKKpipi0::decay
void decay(EvtParticle *p)
Definition EvtDToKKpipi0.cc:159
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
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition EvtPDL.hh:61
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:100
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
m2
double double * m2
Definition qcdloop1.h:83
t
int t()
Definition t.c:1