Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtGoityRoberts.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: EvtGoityRoberts.cc
12//
13// Description: Routine to decay vector-> scalar scalar
14//
15// Modification history:
16//
17// RYD November 24, 1996 Module created
18//
19//------------------------------------------------------------------------
20//
21#include "EvtGenBase/EvtPatches.hh"
22#include <stdlib.h>
23#include "EvtGenBase/EvtParticle.hh"
24#include "EvtGenBase/EvtGenKine.hh"
25#include "EvtGenBase/EvtPDL.hh"
26#include "EvtGenBase/EvtReport.hh"
27#include "EvtGenModels/EvtGoityRoberts.hh"
28#include "EvtGenBase/EvtTensor4C.hh"
29#include "EvtGenBase/EvtDiracSpinor.hh"
30#include <string>
31#include "EvtGenBase/EvtVector4C.hh"
32
33EvtGoityRoberts::~EvtGoityRoberts() {}
34
35std::string EvtGoityRoberts::getName(){
36
37 return "GOITY_ROBERTS";
38
39}
40
41
42EvtDecayBase* EvtGoityRoberts::clone(){
43
44 return new EvtGoityRoberts;
45
46}
47
48void EvtGoityRoberts::init(){
49
50 // check that there are 0 arguments
51 checkNArg(0);
52 checkNDaug(4);
53
54 checkSpinParent(EvtSpinType::SCALAR);
55 checkSpinDaughter(1,EvtSpinType::SCALAR);
56 checkSpinDaughter(2,EvtSpinType::DIRAC);
57 checkSpinDaughter(3,EvtSpinType::NEUTRINO);
58
59}
60
61
62void EvtGoityRoberts::initProbMax() {
63
64 setProbMax( 3000.0);
65}
66
67void EvtGoityRoberts::decay( EvtParticle *p){
68
69 //added by Lange Jan4,2000
70 static EvtId DST0=EvtPDL::getId("D*0");
71 static EvtId DSTB=EvtPDL::getId("anti-D*0");
72 static EvtId DSTP=EvtPDL::getId("D*+");
73 static EvtId DSTM=EvtPDL::getId("D*-");
74 static EvtId D0=EvtPDL::getId("D0");
75 static EvtId D0B=EvtPDL::getId("anti-D0");
76 static EvtId DP=EvtPDL::getId("D+");
77 static EvtId DM=EvtPDL::getId("D-");
78
79
80
81 EvtId meson=getDaug(0);
82
83 if (meson==DST0||meson==DSTP||meson==DSTM||meson==DSTB) {
84 DecayBDstarpilnuGR(p,getDaug(0),getDaug(2),getDaug(3));
85 }
86 else{
87 if (meson==D0||meson==DP||meson==DM||meson==D0B) {
88 DecayBDpilnuGR(p,getDaug(0),getDaug(2),getDaug(3));
89 }
90 else{
91 report(ERROR,"EvtGen") << "Wrong daugther in EvtGoityRoberts!\n";
92 }
93 }
94 return ;
95}
96
97void EvtGoityRoberts::DecayBDstarpilnuGR(EvtParticle *pb,EvtId ndstar,
0ca57c2f 98 EvtId nlep, EvtId /*nnu*/)
da0e9ce3 99{
100
101 pb->initializePhaseSpace(getNDaug(),getDaugs());
102
103 //added by Lange Jan4,2000
104 static EvtId EM=EvtPDL::getId("e-");
105 static EvtId EP=EvtPDL::getId("e+");
106 static EvtId MUM=EvtPDL::getId("mu-");
107 static EvtId MUP=EvtPDL::getId("mu+");
108
109 EvtParticle *dstar, *pion, *lepton, *neutrino;
110
111 // pb->makeDaughters(getNDaug(),getDaugs());
112 dstar=pb->getDaug(0);
113 pion=pb->getDaug(1);
114 lepton=pb->getDaug(2);
115 neutrino=pb->getDaug(3);
116
117 EvtVector4C l1, l2, et0, et1, et2;
118
119 EvtVector4R v,vp,p4_pi;
120 double w;
121
122 v.set(1.0,0.0,0.0,0.0); //4-velocity of B meson
123 vp=(1.0/dstar->getP4().mass())*dstar->getP4(); //4-velocity of D*
124 p4_pi=pion->getP4();
125
126 w=v*vp; //four velocity transfere.
127
128 EvtTensor4C omega;
129
130 double mb=EvtPDL::getMeanMass(pb->getId()); //B mass
131 double md=EvtPDL::getMeanMass(ndstar); //D* mass
132
133 EvtComplex dmb(0.0460,-0.5*0.00001); // B*-B mass splitting ?
134 EvtComplex dmd(0.1421,-0.5*0.00006);
135 // The last two sets of numbers should
136 // be correctly calculated from the
137 // dstar and pion charges.
138 double g = 0.5; // EvtAmplitude proportional to these coupling constants
139 double alpha3 = 0.690; // See table I in G&R's paper
140 double alpha1 = -1.430;
141 double alpha2 = -0.140;
142 double f0 = 0.093; // The pion decay constants set to 93 MeV
143
144 EvtComplex dmt3(0.563,-0.5*0.191); // Mass splitting = dmt - iGamma/2
145 EvtComplex dmt1(0.392,-0.5*1.040);
146 EvtComplex dmt2(0.709,-0.5*0.405);
147
148 double betas=0.285; // magic number for meson wave function ground state
149 double betap=0.280; // magic number for meson wave function state "1"
150 double betad=0.260; // magic number for meson wave function state "2"
151 double betasp=betas*betas+betap*betap;
152 double betasd=betas*betas+betad*betad;
153
154 double lambdabar=0.750; //M(0-,1-) - mQ From Goity&Roberts's code
155
156// Isgur&Wise fct
157 double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
158 double xi1= -1.0*sqrt(2.0/3.0)*(
159 lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
160 exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
161 double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
162 pow((2*betas*betap/(betasp)),2.5)*
163 exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
164 double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
165 pow((2*betas*betad/(betasd)),3.5)*
166 exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
167
168 //report(INFO,"EvtGen") <<"rho's:"<<rho1<<rho2<<endl;
169
170 EvtComplex h1nr,h2nr,h3nr,f1nr,f2nr;
171 EvtComplex f3nr,f4nr,f5nr,f6nr,knr,g1nr,g2nr,g3nr,g4nr,g5nr;
172 EvtComplex h1r,h2r,h3r,f1r,f2r,f3r,f4r,f5r,f6r,kr,g1r,g2r,g3r,g4r,g5r;
173 EvtComplex h1,h2,h3,f1,f2,f3,f4,f5,f6,k,g1,g2,g3,g4,g5;
174
175 // Non-resonance part
176 h1nr = -g*xi*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmb));
177 h2nr = -g*xi/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmb));
178 h3nr = -(g*xi/(f0*md))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
179 -EvtComplex((1.0+w)/(p4_pi*vp),0.0));
180
181 f1nr = -(g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb) -
182 1.0/(EvtComplex(p4_pi*vp,0.0)+dmd));
183 f2nr = f1nr*mb/md;
184 f3nr = EvtComplex(0.0);
185 f4nr = EvtComplex(0.0);
186 f5nr = (g*xi/(2*f0*mb*md))*(EvtComplex(1.0,0.0)
187 +(p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmb));
188 f6nr = (g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
189 -EvtComplex(1.0/(p4_pi*vp),0.0));
190
191 knr = (g*xi/(2*f0))*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmb) +
192 EvtComplex((p4_pi*(v-w*vp))/(p4_pi*vp),0.0));
193
194 g1nr = EvtComplex(0.0);
195 g2nr = EvtComplex(0.0);
196 g3nr = EvtComplex(0.0);
197 g4nr = (g*xi)/(f0*md*EvtComplex(p4_pi*vp));
198 g5nr = EvtComplex(0.0);
199
200 // Resonance part (D** removed by hand - alainb)
201 h1r = -alpha1*rho1*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
202 alpha2*rho2*(p4_pi*(v+2.0*w*v-vp))
203 /(3*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
204 alpha3*xi1*(p4_pi*v)/(f0*mb*md*EvtComplex(p4_pi*v,0.0)+dmt3);
205 h2r = -alpha2*(1+w)*rho2/(3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
206 alpha3*xi1/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
207 h3r = alpha2*rho2*(1+w)/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
208 alpha3*xi1/(f0*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
209
210 f1r = -alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
211 alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
212 f2r = f1r*mb/md;
213 f3r = EvtComplex(0.0);
214 f4r = EvtComplex(0.0);
215 f5r = alpha1*rho1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
216 alpha2*rho2*(p4_pi*(vp-v/3.0-2.0/3.0*w*v))/
217 (2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
218 alpha3*xi1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
219 f6r = alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
220 alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
221
222 kr = -alpha1*rho1*(w-1.0)*(p4_pi*v)/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt1)) -
223 alpha2*rho2*(w-1.0)*(p4_pi*(vp-w*v))
224 /(3*f0*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
225 alpha3*xi1*(p4_pi*(vp-w*v))/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt3));
226
227 g1r = EvtComplex(0.0);
228 g2r = EvtComplex(0.0);
229 g3r = -g2r;
230 g4r = 2.0*alpha2*rho2/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2));
231 g5r = EvtComplex(0.0);
232
233 //Sum
234 h1 = h1nr + h1r;
235 h2 = h2nr + h2r;
236 h3 = h3nr + h3r;
237
238 f1 = f1nr + f1r;
239 f2 = f2nr + f2r;
240 f3 = f3nr + f3r;
241 f4 = f4nr + f4r;
242 f5 = f5nr + f5r;
243 f6 = f6nr + f6r;
244
245 k = knr+kr;
246
247 g1 = g1nr + g1r;
248 g2 = g2nr + g2r;
249 g3 = g3nr + g3r;
250 g4 = g4nr + g4r;
251 g5 = g5nr + g5r;
252
253 EvtTensor4C g_metric;
254 g_metric.setdiag(1.0,-1.0,-1.0,-1.0);
255
256 if (nlep==EM||nlep==MUM){
0ca57c2f 257 omega=EvtComplex(0.0,0.5)*dual(h1*mb*md*EvtGenFunctions::directProd(v,vp)+
258 h2*mb*EvtGenFunctions::directProd(v,p4_pi)+
259 h3*md*EvtGenFunctions::directProd(vp,p4_pi))+
260 f1*mb*EvtGenFunctions::directProd(v,p4_pi)+f2*md*EvtGenFunctions::directProd(vp,p4_pi)+
261 f3*EvtGenFunctions::directProd(p4_pi,p4_pi)+f4*mb*mb*EvtGenFunctions::directProd(v,v)+
262 f5*mb*md*EvtGenFunctions::directProd(vp,v)+f6*mb*EvtGenFunctions::directProd(p4_pi,v)+k*g_metric+
263 EvtComplex(0.0,0.5)*EvtGenFunctions::directProd(dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v),
da0e9ce3 264 (g1*p4_pi+g2*mb*v))+
0ca57c2f 265 EvtComplex(0.0,0.5)*EvtGenFunctions::directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
266 dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v));
da0e9ce3 267
268 l1=EvtLeptonVACurrent(lepton->spParent(0),neutrino->spParentNeutrino());
269 l2=EvtLeptonVACurrent(lepton->spParent(1),neutrino->spParentNeutrino());
270 }
271 else{
272 if (nlep==EP||nlep==MUP){
0ca57c2f 273 omega=EvtComplex(0.0,-0.5)*dual(h1*mb*md*EvtGenFunctions::directProd(v,vp)+
274 h2*mb*EvtGenFunctions::directProd(v,p4_pi)+
275 h3*md*EvtGenFunctions::directProd(vp,p4_pi))+
276 f1*mb*EvtGenFunctions::directProd(v,p4_pi)+f2*md*EvtGenFunctions::directProd(vp,p4_pi)+
277 f3*EvtGenFunctions::directProd(p4_pi,p4_pi)+f4*mb*mb*EvtGenFunctions::directProd(v,v)+
278 f5*mb*md*EvtGenFunctions::directProd(vp,v)+f6*mb*EvtGenFunctions::directProd(p4_pi,v)+k*g_metric+
279 EvtComplex(0.0,-0.5)*EvtGenFunctions::directProd(dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v),
da0e9ce3 280 (g1*p4_pi+g2*mb*v))+
0ca57c2f 281 EvtComplex(0.0,-0.5)*EvtGenFunctions::directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
282 dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v));
da0e9ce3 283
284 l1=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(0));
285 l2=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(1));
286 }
287 else{
288 report(DEBUG,"EvtGen") << "42387dfs878w wrong lepton number\n";
289 }
290 }
291
292 et0=omega.cont2( dstar->epsParent(0).conj() );
293 et1=omega.cont2( dstar->epsParent(1).conj() );
294 et2=omega.cont2( dstar->epsParent(2).conj() );
295
296 vertex(0,0,l1.cont(et0));
297 vertex(0,1,l2.cont(et0));
298
299 vertex(1,0,l1.cont(et1));
300 vertex(1,1,l2.cont(et1));
301
302 vertex(2,0,l1.cont(et2));
303 vertex(2,1,l2.cont(et2));
304
305 return;
306
307}
308
309void EvtGoityRoberts::DecayBDpilnuGR(EvtParticle *pb,EvtId nd,
0ca57c2f 310 EvtId nlep, EvtId /*nnu*/)
da0e9ce3 311
312{
313 //added by Lange Jan4,2000
314 static EvtId EM=EvtPDL::getId("e-");
315 static EvtId EP=EvtPDL::getId("e+");
316 static EvtId MUM=EvtPDL::getId("mu-");
317 static EvtId MUP=EvtPDL::getId("mu+");
318
319 EvtParticle *d, *pion, *lepton, *neutrino;
320
321 pb->initializePhaseSpace(getNDaug(),getDaugs());
322 d=pb->getDaug(0);
323 pion=pb->getDaug(1);
324 lepton=pb->getDaug(2);
325 neutrino=pb->getDaug(3);
326
327 EvtVector4C l1, l2, et0, et1, et2;
328
329 EvtVector4R v,vp,p4_pi;
330 double w;
331
332 v.set(1.0,0.0,0.0,0.0); //4-velocity of B meson
333 vp=(1.0/d->getP4().mass())*d->getP4(); //4-velocity of D
334 p4_pi=pion->getP4(); //4-momentum of pion
335 w=v*vp; //four velocity transfer.
336
337 double mb=EvtPDL::getMeanMass(pb->getId()); //B mass
338 double md=EvtPDL::getMeanMass(nd); //D* mass
339 EvtComplex dmb(0.0460,-0.5*0.00001); //B mass splitting ?
340 //The last two numbers should be
341 //correctly calculated from the
342 //dstar and pion particle number.
343
344 double g = 0.5; // Amplitude proportional to these coupling constants
345 double alpha3 = 0.690; // See table I in G&R's paper
346 double alpha1 = -1.430;
347 double alpha2 = -0.140;
348 double f0=0.093; // The pion decay constant set to 93 MeV
349
350 EvtComplex dmt3(0.563,-0.5*0.191); // Mass splitting = dmt - iGamma/2
351 EvtComplex dmt1(0.392,-0.5*1.040);
352 EvtComplex dmt2(0.709,-0.5*0.405);
353
354 double betas=0.285; // magic number for meson wave function ground state
355 double betap=0.280; // magic number for meson wave function state "1"
356 double betad=0.260; // magic number for meson wave function state "2"
357 double betasp=betas*betas+betap*betap;
358 double betasd=betas*betas+betad*betad;
359
360 double lambdabar=0.750; //M(0-,1-) - mQ From Goity&Roberts's code
361
362 // Isgur&Wise fct
363 double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
364 double xi1= -1.0*sqrt(2.0/3.0)*(lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
365 exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
366 double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
367 pow((2*betas*betap/(betasp)),2.5)*
368 exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
369 double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
370 pow((2*betas*betad/(betasd)),3.5)*
371 exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
372
373 EvtComplex h,a1,a2,a3;
374 EvtComplex hnr,a1nr,a2nr,a3nr;
375 EvtComplex hr,a1r,a2r,a3r;
376
377// Non-resonance part (D* and D** removed by hand - alainb)
378 hnr = g*xi*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb*md);
379 a1nr= -1.0*g*xi*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0);
380 a2nr= g*xi*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb);
381 a3nr=EvtComplex(0.0,0.0);
382
383// Resonance part (D** remove by hand - alainb)
384 hr = alpha2*rho2*(w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0*mb*md) +
385 alpha3*xi1*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb*md);
386 a1r= -1.0*alpha2*rho2*(w*w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0) -
387 alpha3*xi1*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0);
388 a2r= alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*mb) +
389 alpha2*rho2*(0.5*p4_pi*(w*vp-v)+p4_pi*(vp-w*v))/
390 (3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
391 alpha3*xi1*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb);
392 a3r= -1.0*alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*md) -
393 alpha2*rho2*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmt2))/(2*f0*md);
394
395// Sum
396 h=hnr+hr;
397 a1=a1nr+a1r;
398 a2=a2nr+a2r;
399 a3=a3nr+a3r;
400
401 EvtVector4C omega;
402
403 if ( nlep==EM|| nlep==MUM ) {
0ca57c2f 404 omega=EvtComplex(0.0,-1.0)*h*mb*md*dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v)+
da0e9ce3 405 a1*p4_pi+a2*mb*v+a3*md*vp;
406 l1=EvtLeptonVACurrent(
407 lepton->spParent(0),neutrino->spParentNeutrino());
408 l2=EvtLeptonVACurrent(
409 lepton->spParent(1),neutrino->spParentNeutrino());
410 }
411 else{
412 if ( nlep==EP|| nlep==MUP ) {
0ca57c2f 413 omega=EvtComplex(0.0,1.0)*h*mb*md*dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v)+
da0e9ce3 414 a1*p4_pi+a2*mb*v+a3*md*vp;
415 l1=EvtLeptonVACurrent(
416 neutrino->spParentNeutrino(),lepton->spParent(0));
417 l2=EvtLeptonVACurrent(
418 neutrino->spParentNeutrino(),lepton->spParent(1));
419 }
420 else{
421 report(ERROR,"EvtGen") << "42387dfs878w wrong lepton number\n";
422 }
423 }
424
425 vertex(0,l1*omega);
426 vertex(1,l2*omega);
427
428return;
429
430}
431