Line data Source code
1 : //--------------------------------------------------------------------------
2 : //
3 : // Environment:
4 : // This software is part of the EvtGen package developed jointly
5 : // for the BaBar and CLEO collaborations. If you use all or part
6 : // of it, please give an appropriate acknowledgement.
7 : //
8 : // Copyright Information: See EvtGen/COPYRIGHT
9 : // Copyright (C) 1998 Caltech, UCSB
10 : //
11 : // Module: EvtBHadronic.cc
12 : //
13 : // Description: Model for hadronic B decays. Uses naive factorization.
14 : //
15 : // Modification history:
16 : //
17 : // RYD June 14, 1997 Module created
18 : //
19 : //------------------------------------------------------------------------
20 : //
21 : #include "EvtGenBase/EvtPatches.hh"
22 : #include <stdlib.h>
23 : #include "EvtGenBase/EvtParticle.hh"
24 : #include "EvtGenModels/EvtISGW2FF.hh"
25 : #include "EvtGenBase/EvtGenKine.hh"
26 : #include "EvtGenBase/EvtPDL.hh"
27 : #include "EvtGenModels/EvtISGW2FF.hh"
28 : #include "EvtGenBase/EvtTensor4C.hh"
29 : #include "EvtGenBase/EvtVector4C.hh"
30 : #include "EvtGenBase/EvtVector4R.hh"
31 : #include "EvtGenModels/EvtBHadronic.hh"
32 : #include "EvtGenBase/EvtReport.hh"
33 : #include <string>
34 : using std::endl;
35 :
36 0 : EvtBHadronic::~EvtBHadronic() {}
37 :
38 : std::string EvtBHadronic::getName(){
39 :
40 0 : return "BHADRONIC";
41 :
42 : }
43 :
44 :
45 : EvtDecayBase* EvtBHadronic::clone(){
46 :
47 0 : return new EvtBHadronic;
48 :
49 0 : }
50 :
51 :
52 : void EvtBHadronic::init(){
53 :
54 : // check that there are 2 argument
55 0 : checkNArg(2);
56 :
57 :
58 0 : }
59 :
60 : void EvtBHadronic::decay( EvtParticle *p){
61 :
62 : //added by Lange Jan4,2000
63 0 : static EvtId B0=EvtPDL::getId("B0");
64 0 : static EvtId D0=EvtPDL::getId("D0");
65 0 : static EvtId DST0=EvtPDL::getId("D*0");
66 0 : static EvtId D3P10=EvtPDL::getId("D'_10");
67 0 : static EvtId D3P20=EvtPDL::getId("D_2*0");
68 0 : static EvtId D3P00=EvtPDL::getId("D_0*0");
69 0 : static EvtId D1P10=EvtPDL::getId("D_10");
70 :
71 0 : static EvtISGW2FF ffmodel;
72 :
73 0 : p->initializePhaseSpace(getNDaug(),getDaugs());
74 :
75 0 : EvtVector4R p4[MAX_DAUG];
76 0 : p->mass();
77 :
78 :
79 : int i,j;
80 :
81 0 : for ( i=0; i<getNDaug(); i++) {
82 0 : p4[i]=p->getDaug(i)->getP4();
83 : }
84 :
85 : int bcurrent,wcurrent;
86 : int nbcurrent=0;
87 : int nwcurrent=0;
88 :
89 0 : bcurrent=(int)getArg(0);
90 0 : wcurrent=(int)getArg(1);
91 :
92 0 : EvtVector4C jb[5],jw[5];
93 0 : EvtTensor4C g,tds;
94 :
95 0 : EvtVector4R p4b;
96 0 : p4b.set(p->mass(),0.0,0.0,0.0);
97 :
98 0 : EvtVector4R q;
99 : double q2;
100 :
101 0 : EvtComplex ep_meson_bb[5];
102 0 : double f,gf,ap,am;
103 0 : double fp,fm;
104 0 : double hf,kf,bp,bm;
105 0 : EvtVector4R pp,pm;
106 0 : EvtVector4C ep_meson_b[5];
107 :
108 0 : switch (bcurrent) {
109 :
110 : // D0
111 : case 1:
112 0 : q=p4b-p4[0];
113 0 : q2=q*q;
114 : nbcurrent=1;
115 0 : ffmodel.getscalarff(B0,D0,EvtPDL::getMeanMass(D0),
116 0 : EvtPDL::getMeanMass(getDaugs()[1]),&fp,&fm);
117 0 : jb[0]=EvtVector4C(fp*(p4b+p4[0])-fm*q);
118 0 : break;
119 : // D*
120 : case 2:
121 0 : q=p4b-p4[0];
122 0 : q2=q*q;
123 : nbcurrent=3;
124 0 : ffmodel.getvectorff(B0,DST0,EvtPDL::getMeanMass(DST0),q2,&f,&gf,&ap,&am);
125 :
126 0 : g.setdiag(1.0,-1.0,-1.0,-1.0);
127 0 : tds = -f*g
128 0 : -ap*(EvtGenFunctions::directProd(p4b,p4b)+EvtGenFunctions::directProd(p4b,p4[0]))
129 0 : -gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
130 0 : -am*((EvtGenFunctions::directProd(p4b,p4b)-
131 0 : EvtGenFunctions::directProd(p4b,p4[0])));
132 0 : jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
133 0 : jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
134 0 : jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
135 0 : break;
136 : // D1
137 : case 3:
138 0 : q=p4b-p4[0];
139 0 : q2=q*q;
140 : nbcurrent=3;
141 0 : ffmodel.getvectorff(B0,D3P10,EvtPDL::getMeanMass(D3P10),q2,&f,&gf,&ap,&am);
142 :
143 0 : g.setdiag(1.0,-1.0,-1.0,-1.0);
144 0 : tds = -f*g
145 0 : -ap*(EvtGenFunctions::directProd(p4b,p4b)+
146 0 : EvtGenFunctions::directProd(p4b,p4[0]))
147 0 : -gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
148 0 : -am*((EvtGenFunctions::directProd(p4b,p4b)-
149 0 : EvtGenFunctions::directProd(p4b,p4[0])));
150 0 : jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
151 0 : jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
152 0 : jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
153 :
154 0 : break;
155 : // D2*
156 : case 4:
157 0 : q=p4b-p4[0];
158 0 : q2=q*q;
159 : nbcurrent=5;
160 0 : ffmodel.gettensorff(B0,D3P20,EvtPDL::getMeanMass(D3P20),q2,&hf,&kf,&bp,&bm);
161 0 : g.setdiag(1.0,-1.0,-1.0,-1.0);
162 :
163 0 : ep_meson_b[0] = ((p->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
164 0 : ep_meson_b[1] = ((p->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
165 0 : ep_meson_b[2] = ((p->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
166 0 : ep_meson_b[3] = ((p->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
167 0 : ep_meson_b[4] = ((p->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();
168 :
169 0 : pp=p4b+p4[0];
170 0 : pm=p4b-p4[0];
171 :
172 0 : ep_meson_bb[0]=ep_meson_b[0]*(p4b);
173 0 : ep_meson_bb[1]=ep_meson_b[1]*(p4b);
174 0 : ep_meson_bb[2]=ep_meson_b[2]*(p4b);
175 0 : ep_meson_bb[3]=ep_meson_b[3]*(p4b);
176 0 : ep_meson_bb[4]=ep_meson_b[4]*(p4b);
177 :
178 0 : jb[0]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[0])
179 0 : -kf*ep_meson_b[0]
180 0 : -bp*ep_meson_bb[0]*pp-bm*ep_meson_bb[0]*pm;
181 :
182 0 : jb[1]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[1])
183 0 : -kf*ep_meson_b[1]
184 0 : -bp*ep_meson_bb[1]*pp-bm*ep_meson_bb[1]*pm;
185 :
186 0 : jb[2]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[2])
187 0 : -kf*ep_meson_b[2]
188 0 : -bp*ep_meson_bb[2]*pp-bm*ep_meson_bb[2]*pm;
189 :
190 0 : jb[3]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[3])
191 0 : -kf*ep_meson_b[3]
192 0 : -bp*ep_meson_bb[3]*pp-bm*ep_meson_bb[3]*pm;
193 :
194 0 : jb[4]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[4])
195 0 : -kf*ep_meson_b[4]
196 0 : -bp*ep_meson_bb[4]*pp-bm*ep_meson_bb[4]*pm;
197 0 : break;
198 : // D_0*
199 : case 5:
200 0 : q=p4b-p4[0];
201 0 : q2=q*q;
202 0 : double f,gf,ap,am;
203 : nbcurrent=3;
204 0 : ffmodel.getvectorff(B0,D1P10,EvtPDL::getMeanMass(D1P10),q2,&f,&gf,&ap,&am);
205 0 : g.setdiag(1.0,-1.0,-1.0,-1.0);
206 0 : tds = -f*g
207 0 : -ap*(EvtGenFunctions::directProd(p4b,p4b)+
208 0 : EvtGenFunctions::directProd(p4b,p4[0]))
209 0 : +gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
210 0 : -am*((EvtGenFunctions::directProd(p4b,p4b)-
211 0 : EvtGenFunctions::directProd(p4b,p4[0])));
212 0 : jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
213 0 : jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
214 0 : jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
215 0 : break;
216 : // D_1
217 : case 6:
218 0 : q=p4b-p4[0];
219 0 : q2=q*q;
220 : nbcurrent=1;
221 0 : ffmodel.getscalarff(B0,D3P00,EvtPDL::getMeanMass(D3P00),q2,&fp,&fm);
222 0 : jb[0]=fp*(p4b+p4[0])+fm*q;
223 0 : break;
224 : default:
225 0 : report(Severity::Error,"EvtGen")<<"In EvtBHadronic, unknown hadronic current."<<endl;
226 :
227 0 : }
228 :
229 : double norm;
230 :
231 0 : switch (wcurrent) {
232 :
233 :
234 : case 1: case 3: case 4:
235 : nwcurrent=1;
236 0 : jw[0]=p4[getNDaug()-1];
237 0 : break;
238 :
239 : case 2: case 5: case 6:
240 : nwcurrent=3;
241 0 : norm=1.0/sqrt(p4[1].get(0)*p4[1].get(0)/p4[1].mass2()-1.0);
242 0 : jw[0]=norm*p->getDaug(getNDaug()-1)->epsParent(0);
243 0 : jw[1]=norm*p->getDaug(getNDaug()-1)->epsParent(1);
244 0 : jw[2]=norm*p->getDaug(getNDaug()-1)->epsParent(2);
245 0 : break;
246 :
247 :
248 : default:
249 0 : report(Severity::Error,"EvtGen")<<"In EvtBHadronic, unknown W current."<<endl;
250 :
251 : }
252 :
253 0 : if (nbcurrent==1&&nwcurrent==1){
254 0 : vertex(jb[0]*jw[0]);
255 0 : return;
256 : }
257 :
258 0 : if (nbcurrent==1){
259 0 : for(j=0;j<nwcurrent;j++){
260 0 : vertex(j,jb[0]*jw[j]);
261 : }
262 0 : return;
263 : }
264 :
265 0 : if (nwcurrent==1){
266 0 : for(j=0;j<nbcurrent;j++){
267 0 : vertex(j,jb[j]*jw[0]);
268 : }
269 0 : return;
270 : }
271 :
272 0 : for(j=0;j<nbcurrent;j++){
273 0 : for(i=0;i<nwcurrent;i++){
274 0 : vertex(j,i,jb[j]*jw[i]);
275 : }
276 : }
277 0 : return;
278 :
279 0 : }
280 :
281 :
282 :
283 :
284 :
285 :
286 :
|
