]>
Commit | Line | Data |
---|---|---|
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 | ||
33 | EvtGoityRoberts::~EvtGoityRoberts() {} | |
34 | ||
35 | std::string EvtGoityRoberts::getName(){ | |
36 | ||
37 | return "GOITY_ROBERTS"; | |
38 | ||
39 | } | |
40 | ||
41 | ||
42 | EvtDecayBase* EvtGoityRoberts::clone(){ | |
43 | ||
44 | return new EvtGoityRoberts; | |
45 | ||
46 | } | |
47 | ||
48 | void 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 | ||
62 | void EvtGoityRoberts::initProbMax() { | |
63 | ||
64 | setProbMax( 3000.0); | |
65 | } | |
66 | ||
67 | void 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 | ||
97 | void EvtGoityRoberts::DecayBDstarpilnuGR(EvtParticle *pb,EvtId ndstar, | |
98 | EvtId nlep, EvtId nnu) | |
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){ | |
257 | omega=EvtComplex(0.0,0.5)*dual(h1*mb*md*directProd(v,vp)+ | |
258 | h2*mb*directProd(v,p4_pi)+ | |
259 | h3*md*directProd(vp,p4_pi))+ | |
260 | f1*mb*directProd(v,p4_pi)+f2*md*directProd(vp,p4_pi)+ | |
261 | f3*directProd(p4_pi,p4_pi)+f4*mb*mb*directProd(v,v)+ | |
262 | f5*mb*md*directProd(vp,v)+f6*mb*directProd(p4_pi,v)+k*g_metric+ | |
263 | EvtComplex(0.0,0.5)*directProd(dual(directProd(vp,p4_pi)).cont2(v), | |
264 | (g1*p4_pi+g2*mb*v))+ | |
265 | EvtComplex(0.0,0.5)*directProd((g3*mb*v+g4*md*vp+g5*p4_pi), | |
266 | dual(directProd(vp,p4_pi)).cont2(v)); | |
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){ | |
273 | omega=EvtComplex(0.0,-0.5)*dual(h1*mb*md*directProd(v,vp)+ | |
274 | h2*mb*directProd(v,p4_pi)+ | |
275 | h3*md*directProd(vp,p4_pi))+ | |
276 | f1*mb*directProd(v,p4_pi)+f2*md*directProd(vp,p4_pi)+ | |
277 | f3*directProd(p4_pi,p4_pi)+f4*mb*mb*directProd(v,v)+ | |
278 | f5*mb*md*directProd(vp,v)+f6*mb*directProd(p4_pi,v)+k*g_metric+ | |
279 | EvtComplex(0.0,-0.5)*directProd(dual(directProd(vp,p4_pi)).cont2(v), | |
280 | (g1*p4_pi+g2*mb*v))+ | |
281 | EvtComplex(0.0,-0.5)*directProd((g3*mb*v+g4*md*vp+g5*p4_pi), | |
282 | dual(directProd(vp,p4_pi)).cont2(v)); | |
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 | ||
309 | void EvtGoityRoberts::DecayBDpilnuGR(EvtParticle *pb,EvtId nd, | |
310 | EvtId nlep, EvtId nnu) | |
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 ) { | |
404 | omega=EvtComplex(0.0,-1.0)*h*mb*md*dual(directProd(vp,p4_pi)).cont2(v)+ | |
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 ) { | |
413 | omega=EvtComplex(0.0,1.0)*h*mb*md*dual(directProd(vp,p4_pi)).cont2(v)+ | |
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 | ||
428 | return; | |
429 | ||
430 | } | |
431 |