Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtBHadronic.cpp
CommitLineData
da0e9ce3 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>
34using std::endl;
35
36EvtBHadronic::~EvtBHadronic() {}
37
38std::string EvtBHadronic::getName(){
39
40 return "BHADRONIC";
41
42}
43
44
45EvtDecayBase* EvtBHadronic::clone(){
46
47 return new EvtBHadronic;
48
49}
50
51
52void EvtBHadronic::init(){
53
54 // check that there are 2 argument
55 checkNArg(2);
56
57
58}
59
60void EvtBHadronic::decay( EvtParticle *p){
61
62 //added by Lange Jan4,2000
63 static EvtId B0=EvtPDL::getId("B0");
64 static EvtId D0=EvtPDL::getId("D0");
65 static EvtId DST0=EvtPDL::getId("D*0");
66 static EvtId D3P10=EvtPDL::getId("D'_10");
67 static EvtId D3P20=EvtPDL::getId("D_2*0");
68 static EvtId D3P00=EvtPDL::getId("D_0*0");
69 static EvtId D1P10=EvtPDL::getId("D_10");
70
71 static EvtISGW2FF ffmodel;
72
73 p->initializePhaseSpace(getNDaug(),getDaugs());
74
75 EvtVector4R p4[MAX_DAUG];
0ca57c2f 76 p->mass();
da0e9ce3 77
78
79 int i,j;
80
81 for ( i=0; i<getNDaug(); i++) {
82 p4[i]=p->getDaug(i)->getP4();
83 }
84
85 int bcurrent,wcurrent;
86 int nbcurrent=0;
87 int nwcurrent=0;
88
89 bcurrent=(int)getArg(0);
90 wcurrent=(int)getArg(1);
91
92 EvtVector4C jb[5],jw[5];
93 EvtTensor4C g,tds;
94
95 EvtVector4R p4b;
96 p4b.set(p->mass(),0.0,0.0,0.0);
97
98 EvtVector4R q;
99 double q2;
100
101 EvtComplex ep_meson_bb[5];
102 double f,gf,ap,am;
103 double fp,fm;
104 double hf,kf,bp,bm;
105 EvtVector4R pp,pm;
106 EvtVector4C ep_meson_b[5];
107
108 switch (bcurrent) {
109
110 // D0
111 case 1:
112 q=p4b-p4[0];
113 q2=q*q;
114 nbcurrent=1;
115 ffmodel.getscalarff(B0,D0,EvtPDL::getMeanMass(D0),
116 EvtPDL::getMeanMass(getDaugs()[1]),&fp,&fm);
117 jb[0]=EvtVector4C(fp*(p4b+p4[0])-fm*q);
118 break;
119 // D*
120 case 2:
121 q=p4b-p4[0];
122 q2=q*q;
123 nbcurrent=3;
124 ffmodel.getvectorff(B0,DST0,EvtPDL::getMeanMass(DST0),q2,&f,&gf,&ap,&am);
125
126 g.setdiag(1.0,-1.0,-1.0,-1.0);
127 tds = -f*g
0ca57c2f 128 -ap*(EvtGenFunctions::directProd(p4b,p4b)+EvtGenFunctions::directProd(p4b,p4[0]))
129 -gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
130 -am*((EvtGenFunctions::directProd(p4b,p4b)-
131 EvtGenFunctions::directProd(p4b,p4[0])));
da0e9ce3 132 jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
133 jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
134 jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
135 break;
136 // D1
137 case 3:
138 q=p4b-p4[0];
139 q2=q*q;
140 nbcurrent=3;
141 ffmodel.getvectorff(B0,D3P10,EvtPDL::getMeanMass(D3P10),q2,&f,&gf,&ap,&am);
142
143 g.setdiag(1.0,-1.0,-1.0,-1.0);
144 tds = -f*g
0ca57c2f 145 -ap*(EvtGenFunctions::directProd(p4b,p4b)+
146 EvtGenFunctions::directProd(p4b,p4[0]))
147 -gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
148 -am*((EvtGenFunctions::directProd(p4b,p4b)-
149 EvtGenFunctions::directProd(p4b,p4[0])));
da0e9ce3 150 jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
151 jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
152 jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
153
154 break;
155 // D2*
156 case 4:
157 q=p4b-p4[0];
158 q2=q*q;
159 nbcurrent=5;
160 ffmodel.gettensorff(B0,D3P20,EvtPDL::getMeanMass(D3P20),q2,&hf,&kf,&bp,&bm);
161 g.setdiag(1.0,-1.0,-1.0,-1.0);
162
163 ep_meson_b[0] = ((p->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
164 ep_meson_b[1] = ((p->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
165 ep_meson_b[2] = ((p->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
166 ep_meson_b[3] = ((p->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
167 ep_meson_b[4] = ((p->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();
168
169 pp=p4b+p4[0];
170 pm=p4b-p4[0];
171
172 ep_meson_bb[0]=ep_meson_b[0]*(p4b);
173 ep_meson_bb[1]=ep_meson_b[1]*(p4b);
174 ep_meson_bb[2]=ep_meson_b[2]*(p4b);
175 ep_meson_bb[3]=ep_meson_b[3]*(p4b);
176 ep_meson_bb[4]=ep_meson_b[4]*(p4b);
177
0ca57c2f 178 jb[0]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[0])
da0e9ce3 179 -kf*ep_meson_b[0]
180 -bp*ep_meson_bb[0]*pp-bm*ep_meson_bb[0]*pm;
181
0ca57c2f 182 jb[1]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[1])
da0e9ce3 183 -kf*ep_meson_b[1]
184 -bp*ep_meson_bb[1]*pp-bm*ep_meson_bb[1]*pm;
185
0ca57c2f 186 jb[2]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[2])
da0e9ce3 187 -kf*ep_meson_b[2]
188 -bp*ep_meson_bb[2]*pp-bm*ep_meson_bb[2]*pm;
189
0ca57c2f 190 jb[3]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[3])
da0e9ce3 191 -kf*ep_meson_b[3]
192 -bp*ep_meson_bb[3]*pp-bm*ep_meson_bb[3]*pm;
193
0ca57c2f 194 jb[4]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[4])
da0e9ce3 195 -kf*ep_meson_b[4]
196 -bp*ep_meson_bb[4]*pp-bm*ep_meson_bb[4]*pm;
197 break;
198 // D_0*
199 case 5:
200 q=p4b-p4[0];
201 q2=q*q;
202 double f,gf,ap,am;
203 nbcurrent=3;
204 ffmodel.getvectorff(B0,D1P10,EvtPDL::getMeanMass(D1P10),q2,&f,&gf,&ap,&am);
205 g.setdiag(1.0,-1.0,-1.0,-1.0);
206 tds = -f*g
0ca57c2f 207 -ap*(EvtGenFunctions::directProd(p4b,p4b)+
208 EvtGenFunctions::directProd(p4b,p4[0]))
209 +gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
210 -am*((EvtGenFunctions::directProd(p4b,p4b)-
211 EvtGenFunctions::directProd(p4b,p4[0])));
da0e9ce3 212 jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
213 jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
214 jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
215 break;
216 // D_1
217 case 6:
218 q=p4b-p4[0];
219 q2=q*q;
220 nbcurrent=1;
221 ffmodel.getscalarff(B0,D3P00,EvtPDL::getMeanMass(D3P00),q2,&fp,&fm);
222 jb[0]=fp*(p4b+p4[0])+fm*q;
223 break;
224 default:
225 report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown hadronic current."<<endl;
226
227 }
228
229 double norm;
230
231 switch (wcurrent) {
232
233
234 case 1: case 3: case 4:
235 nwcurrent=1;
236 jw[0]=p4[getNDaug()-1];
237 break;
238
239 case 2: case 5: case 6:
240 nwcurrent=3;
241 norm=1.0/sqrt(p4[1].get(0)*p4[1].get(0)/p4[1].mass2()-1.0);
242 jw[0]=norm*p->getDaug(getNDaug()-1)->epsParent(0);
243 jw[1]=norm*p->getDaug(getNDaug()-1)->epsParent(1);
244 jw[2]=norm*p->getDaug(getNDaug()-1)->epsParent(2);
245 break;
246
247
248 default:
249 report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown W current."<<endl;
250
251 }
252
253 if (nbcurrent==1&&nwcurrent==1){
254 vertex(jb[0]*jw[0]);
255 return;
256 }
257
258 if (nbcurrent==1){
259 for(j=0;j<nwcurrent;j++){
260 vertex(j,jb[0]*jw[j]);
261 }
262 return;
263 }
264
265 if (nwcurrent==1){
266 for(j=0;j<nbcurrent;j++){
267 vertex(j,jb[j]*jw[0]);
268 }
269 return;
270 }
271
272 for(j=0;j<nbcurrent;j++){
273 for(i=0;i<nwcurrent;i++){
274 vertex(j,i,jb[j]*jw[i]);
275 }
276 }
277 return;
278
279}
280
281
282
283
284
285
286