BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToPiPi0Etap.cc
Go to the documentation of this file.
1//------------------------------------------------------------------------
2// Environment:
3// This software is part of the EvtGen package developed jointly
4// for the BaBar and CLEO collaborations. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/COPYRIGHT
8// Copyright (C) 1998 Caltech, UCSB
9//
10// Module: EvtD0ToEtapipi.cc
11//
12// Description: Routine to handle three-body decays of D0/D0_bar or D+/D-
13//
14// Modification history:
15//
16// Liaoyuan Dong Jul 3 23:56:42 2023 Module created
17// Xian Kui Oct 5 20:45:00 2023 DToPiPi0Etap
18//------------------------------------------------------------------------
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 "EvtGenBase/EvtVector4R.hh"
26#include "EvtGenBase/EvtComplex.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
28#include "EvtGenModels/EvtDToPiPi0Etap.hh"
29#include <stdlib.h>
30using namespace std;
31
32EvtDToPiPi0Etap::~EvtDToPiPi0Etap() {}
33
34void EvtDToPiPi0Etap::getName(std::string& model_name){
35 model_name = "DToPiPi0Etap";
36}
37
41
42void EvtDToPiPi0Etap::init(){
43
44 checkNArg(0);
45 checkNDaug(3);
46 checkSpinParent(EvtSpinType::SCALAR);
47 checkSpinDaughter(0,EvtSpinType::SCALAR);
48 checkSpinDaughter(1,EvtSpinType::SCALAR);
49 checkSpinDaughter(2,EvtSpinType::SCALAR);
50
51 modetype[0] = 12;
52 modetype[1] = 12;
53 rho[0] = 1.0; // rho+
54 phi[0] = 0.0;
55 rho[1] = 8.3996e+00; // (Pi+Pi0) D-wave
56 phi[1] = 4.8088e+00;
57
58 mass[0] = 0.77511;
59 width[0] = 0.1491;
60
61 mass[1] = 0.27454719;
62 width[1] = 0.1491; // TBD
63
64 massDp = 1.86966;
65 massD0 = 1.86484;
66
67 massEtap = 0.95778;
68 massPip = 0.13957;
69 massPi0 = 0.13498;
70
71 mathPi = 3.1415926;
72
73 // GS
74 GS1 = 0.6366197830;
75 GS2 = 0.0186018247;
76 GS3 = 0.1591549431; // 1/(2*mathPi)
77 GS4 = 0.0062006081; // massPip*massPip/mathPi
78
79 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
80
81 for (int i=0; i<4; i++) {
82 for (int j=0; j<4; j++) {
83 G[i][j] = GG[i][j];
84 }
85 }
86}
87
88void EvtDToPiPi0Etap::initProbMax() {
89 setProbMax(48.9); // The Max Prob = 48.895731256754906724
90}
91
92void EvtDToPiPi0Etap::decay(EvtParticle *p){
93
94 // This piece of code could in principle be used to calculate maximum
95 // probablity on fly. But as it uses high number of points and model
96 // deals with single final state, we keep hardcoded number for now rather
97 // than adapting code to work here.
98
99/*
100 double maxprob = 0.0;
101 for(int ir=0;ir<=60000000;ir++){
102 p->initializePhaseSpace(getNDaug(),getDaugs());
103 EvtVector4R D1 = p->getDaug(0)->getP4();
104 EvtVector4R D2 = p->getDaug(1)->getP4();
105 EvtVector4R D3 = p->getDaug(2)->getP4();
106
107 double P1[4], P2[4], P3[4];
108 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
109 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
110 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
111
112 double value;
113 int g0[2]={1,0};
114 int spin[2]={1,2};
115 int nstates=2;
116
117 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
118 if(value>maxprob) {
119 maxprob=value;
120 cout << "ir = " << ir << " maxProb= " << value << endl;
121 }
122 }
123 cout << "The maxProb = " << maxprob << endl;
124*/
125
126 p->initializePhaseSpace(getNDaug(),getDaugs());
127 EvtVector4R D1 = p->getDaug(0)->getP4();
128 EvtVector4R D2 = p->getDaug(1)->getP4();
129 EvtVector4R D3 = p->getDaug(2)->getP4();
130
131 double P1[4], P2[4], P3[4];
132
133 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
134 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
135 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
136
137 double value;
138 int g0[2]={1,0};
139 int spin[2]={1,2};
140 int nstates=2;
141 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
142
143 setProb(value);
144 return ;
145}
146
147/* *
148 * Base
149 * */
150
151void EvtDToPiPi0Etap::Com_Multi(double a1[2], double a2[2], double res[2])
152{
153 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
154 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
155}
156
157void EvtDToPiPi0Etap::Com_Divide(double a1[2], double a2[2], double res[2])
158{
159 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
160 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
161 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
162}
163
164double EvtDToPiPi0Etap::SCADot(double a1[4], double a2[4])
165{
166 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
167 return _cal;
168}
169
170double EvtDToPiPi0Etap::barrier(int l, double sa, double sb, double sc, double r, double mass)
171{
172 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
173 if(q < 0) q = -q;
174 double z;
175 z = q*r*r;
176 double sa0;
177 sa0 = mass*mass;
178 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
179 if(q0 < 0) q0 = -1*q0;
180 double z0 = q0*r*r;
181 double F = 0.0;
182 if(l == 0) F = 1;
183 if(l == 1) F = sqrt((1+z0)/(1+z));
184 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
185 return F;
186}
187
188/* *
189 * Spin
190 * */
191
192void EvtDToPiPi0Etap::calt1(double daug1[4], double daug2[4], double t1[4])
193{
194 double p, pq, tmp;
195 double pa[4], qa[4];
196 for(int i=0; i<4; i++) {
197 pa[i] = daug1[i] + daug2[i];
198 qa[i] = daug1[i] - daug2[i];
199 }
200 p = SCADot(pa,pa);
201 pq = SCADot(pa,qa);
202 tmp = pq/p;
203 for(int i=0; i<4; i++) {
204 t1[i] = qa[i] - tmp*pa[i];
205 }
206}
207
208void EvtDToPiPi0Etap::calt2(double daug1[4], double daug2[4], double t2[4][4])
209{
210 double p, r;
211 double pa[4], t1[4];
212 calt1(daug1,daug2,t1);
213 r = SCADot(t1,t1)/3.0;
214 for(int i=0; i<4; i++) {
215 pa[i] = daug1[i] + daug2[i];
216 }
217 p = SCADot(pa,pa);
218 for(int i=0; i<4; i++) {
219 for(int j=0; j<4; j++) {
220 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
221 }
222 }
223}
224
225/* *
226 * Propagator
227 * */
228
229double EvtDToPiPi0Etap::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
230{
231 double widm = 0.;
232 double m = sqrt(sa);
233 double tmp = sb-sc;
234 double tmp1 = sa+tmp;
235 double q = 0.25*tmp1*tmp1/sa-sb;
236 if(q<0) q = -q;
237 double tmp2 = mass2+tmp;
238 double q0 = 0.25*tmp2*tmp2/mass2-sb;
239 if(q0<0) q0 = -q0;
240 double z = q*r2;
241 double z0 = q0*r2;
242 double t = q/q0;
243 if(l == 0) {widm = sqrt(t)*mass/m;}
244 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
245 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
246 return widm;
247}
248
249double EvtDToPiPi0Etap::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
250{
251 double widm = 0.;
252 double m = sqrt(sa);
253 double tmp = sb-sc;
254 double tmp1 = sa+tmp;
255 double q = 0.25*tmp1*tmp1/sa-sb;
256 if(q<0) q = -q;
257 double tmp2 = mass2+tmp;
258 double q0 = 0.25*tmp2*tmp2/mass2-sb;
259 if(q0<0) q0 = -q0;
260 double z = q*r2;
261 double z0 = q0*r2;
262 double F = (1+z0)/(1+z);
263 double t = q/q0;
264 widm = t*sqrt(t)*mass/m*F;
265 return widm;
266}
267
268void EvtDToPiPi0Etap::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
269{
270 double a[2], b[2];
271 a[0] = 1;
272 a[1] = 0;
273
274 b[0] = mass2-sa;
275 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
276 Com_Divide(a,b,prop);
277}
278
279void EvtDToPiPi0Etap::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
280{
281 double a[2], b[2];
282 double tmp = sb-sc;
283 double tmp1 = sa+tmp;
284 double q2 = 0.25*tmp1*tmp1/sa-sb;
285 if(q2<0) q2 = 1e-16;
286
287 double tmp2 = mass2+tmp;
288 double q02 = 0.25*tmp2*tmp2/mass2-sb;
289 if(q02<0) q02 = 1e-16;
290
291 double q = sqrt(q2);
292 double q0 = sqrt(q02);
293 double m = sqrt(sa);
294 double q03 = q0*q02;
295 double tmp3 = log(mass+2*q0) + 1.292621885; // ln(massPip + massPi0) = -1.292621885;
296
297 double h = GS1*q/m*(log(m+2*q) + 1.292621885);
298 double h0 = GS1*q0/mass*tmp3;
299 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
300 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
301 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
302
303 a[0] = 1.0+d*width/mass;
304 a[1] = 0.0;
305 b[0] = mass2-sa+width*f;
306 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
307 Com_Divide(a,b,prop);
308}
309
310double EvtDToPiPi0Etap::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass)
311{
312 double pR[4], pD[4];
313 double temp_PDF, v_re;
314 temp_PDF = 0.0;
315 v_re = 0.0;
316 double B[2], s1, s2, s3, sR, sD;
317 for(int i=0; i<4; i++){
318 pR[i] = P1[i] + P2[i];
319 pD[i] = pR[i] + P3[i];
320 }
321 s1 = SCADot(P1,P1);
322 s2 = SCADot(P2,P2);
323 s3 = SCADot(P3,P3);
324 sR = SCADot(pR,pR);
325 sD = SCADot(pD,pD);
326 int G[4][4];
327 for(int i=0; i!=4; i++){
328 for(int j=0; j!=4; j++){
329 if(i==j){
330 if(i==0) G[i][j] = 1;
331 else G[i][j] = -1;
332 }
333 else G[i][j] = 0;
334 }
335 }
336
337 if(Ang == 0){
338 B[0] = 1;
339 B[1] = 1;
340 temp_PDF = 1;
341 }
342
343 if(Ang == 1){
344 B[0] = barrier(1,sR,s1,s2,3.0,mass);
345 B[1] = barrier(1,sD,sR,s3,5.0,massDp);
346 double t1[4], T1[4];
347 calt1(P1,P2,t1);
348 calt1(pR,P3,T1);
349 temp_PDF = 0;
350 for(int i=0; i<4; i++){
351 temp_PDF += t1[i]*T1[i]*G[i][i];
352 }
353 }
354
355 if(Ang == 2){
356 B[0] = barrier(2,sR,s1,s2,3.0,mass);
357 B[1] = barrier(2,sD,sR,s3,5.0,massDp);
358 double t2[4][4], T2[4][4];
359 calt2(P1,P2,t2);
360 calt2(pR,P3,T2);
361 temp_PDF = 0;
362 for(int i=0; i<4; i++){
363 for(int j=0; j<4; j++){
364 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
365 }
366 }
367 }
368
369 v_re = temp_PDF*B[0]*B[1];
370 return v_re;
371}
372
373/* *
374 * Calculation
375 * */
376
377void EvtDToPiPi0Etap::calEva(
378 double* Pip, double* Pi0, double* Etap, double *mass, double *width,
379 double *amp, double *phase, int* g0, int* spin, int* modetype, int nstates, double & Result)
380{
381 double P12[4], P23[4], P13[4];
382 double coEff[2], tempPDF[2], ultPDF[2];
383 double rRes = 9;
384
385 double sPi0, sEtap, sPip;
386 double s12, s13, s23;
387
388 for(int i=0; i<4; i++){
389 P12[i] = Pip[i] + Pi0[i];
390 P13[i] = Pip[i] + Etap[i];
391 P23[i] = Pi0[i] + Etap[i];
392 }
393
394 sPip = SCADot(Pip,Pip);
395 sPi0 = SCADot(Pi0,Pi0);
396 sEtap = SCADot(Etap,Etap);
397
398 s12 = SCADot(P12,P12);
399 s13 = SCADot(P13,P13);
400 s23 = SCADot(P23,P23);
401
402 double prop[2], partAmp, tempAmp[2];
403 double massSq;
404
405 ultPDF[0] = 0;
406 ultPDF[1] = 0;
407
408 tempPDF[0] = 0;
409 tempPDF[1] = 0;
410
411 for(int i = 0; i<nstates; i++) {
412 tempAmp[0] = 0;
413 tempAmp[1] = 0;
414 partAmp = 0;
415
416 massSq = mass[i]*mass[i];
417 coEff[0] = amp[i]*cos(phase[i]);
418 coEff[1] = amp[i]*sin(phase[i]);
419
420 if(modetype[i] == 12){
421 partAmp = DDalitz(Pip, Pi0, Etap, spin[i], mass[i]);
422 if(g0[i]==1) propagatorGS(massSq, mass[i], width[i], s12, sPip, sPi0, rRes, prop);
423 if(g0[i]==0){
424 prop[0] = 1;
425 prop[1] = 0;
426 }
427 tempAmp[0] = partAmp*prop[0];
428 tempAmp[1] = partAmp*prop[1];
429 }
430
431 if(modetype[i] == 13){
432 partAmp = DDalitz(Pip, Etap, Pi0, spin[i], mass[i]);
433 if(g0[i]==1) propagatorRBW(massSq, mass[i], width[i], s13, sPip, sEtap, rRes, spin[i], prop);
434 if(g0[i]==0){
435 prop[0] = 1;
436 prop[1] = 0;
437 }
438 tempAmp[0] = partAmp*prop[0];
439 tempAmp[1] = partAmp*prop[1];
440 }
441
442 if(modetype[i] == 23){
443 partAmp = DDalitz(Pi0, Etap, Pip, spin[i], mass[i]);
444 if(g0[i]==1) propagatorRBW(massSq, mass[i], width[i], s23, sPi0, sEtap, rRes, spin[i], prop);
445 if(g0[i]==0){
446 prop[0] = 1;
447 prop[1] = 0;
448 }
449 tempAmp[0] = partAmp*prop[0];
450 tempAmp[1] = partAmp*prop[1];
451 }
452
453 if(modetype[i] == 1111) {
454 tempAmp[0] = 1;
455 tempAmp[1] = 0;
456 }
457
458 Com_Multi(tempAmp, coEff, tempPDF);
459
460 ultPDF[0] += tempPDF[0];
461 ultPDF[1] += tempPDF[1];
462 }
463
464 double value = ultPDF[0]*ultPDF[0] + ultPDF[1]*ultPDF[1];
465
466 if(value <=0) value = 1e-20;
467 Result = value;
468}
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")
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
t
TTree * t
Definition binning.cxx:23
EvtDToPiPi0Etap::~EvtDToPiPi0Etap
virtual ~EvtDToPiPi0Etap()
Definition EvtDToPiPi0Etap.cc:32
EvtDToPiPi0Etap::decay
void decay(EvtParticle *p)
Definition EvtDToPiPi0Etap.cc:92
EvtDToPiPi0Etap::clone
EvtDecayBase * clone()
Definition EvtDToPiPi0Etap.cc:38
EvtDToPiPi0Etap::getName
void getName(std::string &name)
Definition EvtDToPiPi0Etap.cc:34
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:533
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