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