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: 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 0 : EvtGoityRoberts::~EvtGoityRoberts() {}
34 :
35 : std::string EvtGoityRoberts::getName(){
36 :
37 0 : return "GOITY_ROBERTS";
38 :
39 : }
40 :
41 :
42 : EvtDecayBase* EvtGoityRoberts::clone(){
43 :
44 0 : return new EvtGoityRoberts;
45 :
46 0 : }
47 :
48 : void EvtGoityRoberts::init(){
49 :
50 : // check that there are 0 arguments
51 0 : checkNArg(0);
52 0 : checkNDaug(4);
53 :
54 0 : checkSpinParent(EvtSpinType::SCALAR);
55 0 : checkSpinDaughter(1,EvtSpinType::SCALAR);
56 0 : checkSpinDaughter(2,EvtSpinType::DIRAC);
57 0 : checkSpinDaughter(3,EvtSpinType::NEUTRINO);
58 :
59 0 : }
60 :
61 :
62 : void EvtGoityRoberts::initProbMax() {
63 :
64 0 : setProbMax( 3000.0);
65 0 : }
66 :
67 : void EvtGoityRoberts::decay( EvtParticle *p){
68 :
69 : //added by Lange Jan4,2000
70 0 : static EvtId DST0=EvtPDL::getId("D*0");
71 0 : static EvtId DSTB=EvtPDL::getId("anti-D*0");
72 0 : static EvtId DSTP=EvtPDL::getId("D*+");
73 0 : static EvtId DSTM=EvtPDL::getId("D*-");
74 0 : static EvtId D0=EvtPDL::getId("D0");
75 0 : static EvtId D0B=EvtPDL::getId("anti-D0");
76 0 : static EvtId DP=EvtPDL::getId("D+");
77 0 : static EvtId DM=EvtPDL::getId("D-");
78 :
79 :
80 :
81 0 : EvtId meson=getDaug(0);
82 :
83 0 : if (meson==DST0||meson==DSTP||meson==DSTM||meson==DSTB) {
84 0 : DecayBDstarpilnuGR(p,getDaug(0),getDaug(2),getDaug(3));
85 0 : }
86 : else{
87 0 : if (meson==D0||meson==DP||meson==DM||meson==D0B) {
88 0 : DecayBDpilnuGR(p,getDaug(0),getDaug(2),getDaug(3));
89 0 : }
90 : else{
91 0 : report(Severity::Error,"EvtGen") << "Wrong daugther in EvtGoityRoberts!\n";
92 : }
93 : }
94 : return ;
95 0 : }
96 :
97 : void EvtGoityRoberts::DecayBDstarpilnuGR(EvtParticle *pb,EvtId ndstar,
98 : EvtId nlep, EvtId /*nnu*/)
99 : {
100 :
101 0 : pb->initializePhaseSpace(getNDaug(),getDaugs());
102 :
103 : //added by Lange Jan4,2000
104 0 : static EvtId EM=EvtPDL::getId("e-");
105 0 : static EvtId EP=EvtPDL::getId("e+");
106 0 : static EvtId MUM=EvtPDL::getId("mu-");
107 0 : static EvtId MUP=EvtPDL::getId("mu+");
108 :
109 : EvtParticle *dstar, *pion, *lepton, *neutrino;
110 :
111 : // pb->makeDaughters(getNDaug(),getDaugs());
112 0 : dstar=pb->getDaug(0);
113 0 : pion=pb->getDaug(1);
114 0 : lepton=pb->getDaug(2);
115 0 : neutrino=pb->getDaug(3);
116 :
117 0 : EvtVector4C l1, l2, et0, et1, et2;
118 :
119 0 : EvtVector4R v,vp,p4_pi;
120 : double w;
121 :
122 0 : v.set(1.0,0.0,0.0,0.0); //4-velocity of B meson
123 0 : vp=(1.0/dstar->getP4().mass())*dstar->getP4(); //4-velocity of D*
124 0 : p4_pi=pion->getP4();
125 :
126 0 : w=v*vp; //four velocity transfere.
127 :
128 0 : EvtTensor4C omega;
129 :
130 0 : double mb=EvtPDL::getMeanMass(pb->getId()); //B mass
131 0 : double md=EvtPDL::getMeanMass(ndstar); //D* mass
132 :
133 0 : EvtComplex dmb(0.0460,-0.5*0.00001); // B*-B mass splitting ?
134 0 : 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 0 : EvtComplex dmt3(0.563,-0.5*0.191); // Mass splitting = dmt - iGamma/2
145 0 : EvtComplex dmt1(0.392,-0.5*1.040);
146 0 : 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 0 : double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
158 0 : double xi1= -1.0*sqrt(2.0/3.0)*(
159 0 : lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
160 : exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
161 0 : double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
162 0 : pow((2*betas*betap/(betasp)),2.5)*
163 0 : exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
164 0 : double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
165 0 : pow((2*betas*betad/(betasd)),3.5)*
166 0 : exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
167 :
168 : //report(Severity::Info,"EvtGen") <<"rho's:"<<rho1<<rho2<<endl;
169 :
170 0 : EvtComplex h1nr,h2nr,h3nr,f1nr,f2nr;
171 0 : EvtComplex f3nr,f4nr,f5nr,f6nr,knr,g1nr,g2nr,g3nr,g4nr,g5nr;
172 0 : EvtComplex h1r,h2r,h3r,f1r,f2r,f3r,f4r,f5r,f6r,kr,g1r,g2r,g3r,g4r,g5r;
173 0 : EvtComplex h1,h2,h3,f1,f2,f3,f4,f5,f6,k,g1,g2,g3,g4,g5;
174 :
175 : // Non-resonance part
176 0 : h1nr = -g*xi*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmb));
177 0 : h2nr = -g*xi/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmb));
178 0 : h3nr = -(g*xi/(f0*md))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
179 0 : -EvtComplex((1.0+w)/(p4_pi*vp),0.0));
180 :
181 0 : f1nr = -(g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb) -
182 0 : 1.0/(EvtComplex(p4_pi*vp,0.0)+dmd));
183 0 : f2nr = f1nr*mb/md;
184 0 : f3nr = EvtComplex(0.0);
185 0 : f4nr = EvtComplex(0.0);
186 0 : f5nr = (g*xi/(2*f0*mb*md))*(EvtComplex(1.0,0.0)
187 0 : +(p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmb));
188 0 : f6nr = (g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
189 0 : -EvtComplex(1.0/(p4_pi*vp),0.0));
190 :
191 0 : knr = (g*xi/(2*f0))*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmb) +
192 0 : EvtComplex((p4_pi*(v-w*vp))/(p4_pi*vp),0.0));
193 :
194 0 : g1nr = EvtComplex(0.0);
195 0 : g2nr = EvtComplex(0.0);
196 0 : g3nr = EvtComplex(0.0);
197 0 : g4nr = (g*xi)/(f0*md*EvtComplex(p4_pi*vp));
198 0 : g5nr = EvtComplex(0.0);
199 :
200 : // Resonance part (D** removed by hand - alainb)
201 0 : h1r = -alpha1*rho1*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
202 0 : alpha2*rho2*(p4_pi*(v+2.0*w*v-vp))
203 0 : /(3*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
204 0 : alpha3*xi1*(p4_pi*v)/(f0*mb*md*EvtComplex(p4_pi*v,0.0)+dmt3);
205 0 : h2r = -alpha2*(1+w)*rho2/(3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
206 0 : alpha3*xi1/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
207 0 : h3r = alpha2*rho2*(1+w)/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
208 0 : alpha3*xi1/(f0*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
209 :
210 0 : f1r = -alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
211 0 : alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
212 0 : f2r = f1r*mb/md;
213 0 : f3r = EvtComplex(0.0);
214 0 : f4r = EvtComplex(0.0);
215 0 : f5r = alpha1*rho1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
216 0 : alpha2*rho2*(p4_pi*(vp-v/3.0-2.0/3.0*w*v))/
217 0 : (2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
218 0 : alpha3*xi1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
219 0 : f6r = alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
220 0 : alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
221 :
222 0 : kr = -alpha1*rho1*(w-1.0)*(p4_pi*v)/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt1)) -
223 0 : alpha2*rho2*(w-1.0)*(p4_pi*(vp-w*v))
224 0 : /(3*f0*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
225 0 : alpha3*xi1*(p4_pi*(vp-w*v))/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt3));
226 :
227 0 : g1r = EvtComplex(0.0);
228 0 : g2r = EvtComplex(0.0);
229 0 : g3r = -g2r;
230 0 : g4r = 2.0*alpha2*rho2/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2));
231 0 : g5r = EvtComplex(0.0);
232 :
233 : //Sum
234 0 : h1 = h1nr + h1r;
235 0 : h2 = h2nr + h2r;
236 0 : h3 = h3nr + h3r;
237 :
238 0 : f1 = f1nr + f1r;
239 0 : f2 = f2nr + f2r;
240 0 : f3 = f3nr + f3r;
241 0 : f4 = f4nr + f4r;
242 0 : f5 = f5nr + f5r;
243 0 : f6 = f6nr + f6r;
244 :
245 0 : k = knr+kr;
246 :
247 0 : g1 = g1nr + g1r;
248 0 : g2 = g2nr + g2r;
249 0 : g3 = g3nr + g3r;
250 0 : g4 = g4nr + g4r;
251 0 : g5 = g5nr + g5r;
252 :
253 0 : EvtTensor4C g_metric;
254 0 : g_metric.setdiag(1.0,-1.0,-1.0,-1.0);
255 :
256 0 : if (nlep==EM||nlep==MUM){
257 0 : omega=EvtComplex(0.0,0.5)*dual(h1*mb*md*EvtGenFunctions::directProd(v,vp)+
258 0 : h2*mb*EvtGenFunctions::directProd(v,p4_pi)+
259 0 : h3*md*EvtGenFunctions::directProd(vp,p4_pi))+
260 0 : f1*mb*EvtGenFunctions::directProd(v,p4_pi)+f2*md*EvtGenFunctions::directProd(vp,p4_pi)+
261 0 : f3*EvtGenFunctions::directProd(p4_pi,p4_pi)+f4*mb*mb*EvtGenFunctions::directProd(v,v)+
262 0 : f5*mb*md*EvtGenFunctions::directProd(vp,v)+f6*mb*EvtGenFunctions::directProd(p4_pi,v)+k*g_metric+
263 0 : EvtComplex(0.0,0.5)*EvtGenFunctions::directProd(dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v),
264 0 : (g1*p4_pi+g2*mb*v))+
265 0 : EvtComplex(0.0,0.5)*EvtGenFunctions::directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
266 0 : dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v));
267 :
268 0 : l1=EvtLeptonVACurrent(lepton->spParent(0),neutrino->spParentNeutrino());
269 0 : l2=EvtLeptonVACurrent(lepton->spParent(1),neutrino->spParentNeutrino());
270 0 : }
271 : else{
272 0 : if (nlep==EP||nlep==MUP){
273 0 : omega=EvtComplex(0.0,-0.5)*dual(h1*mb*md*EvtGenFunctions::directProd(v,vp)+
274 0 : h2*mb*EvtGenFunctions::directProd(v,p4_pi)+
275 0 : h3*md*EvtGenFunctions::directProd(vp,p4_pi))+
276 0 : f1*mb*EvtGenFunctions::directProd(v,p4_pi)+f2*md*EvtGenFunctions::directProd(vp,p4_pi)+
277 0 : f3*EvtGenFunctions::directProd(p4_pi,p4_pi)+f4*mb*mb*EvtGenFunctions::directProd(v,v)+
278 0 : f5*mb*md*EvtGenFunctions::directProd(vp,v)+f6*mb*EvtGenFunctions::directProd(p4_pi,v)+k*g_metric+
279 0 : EvtComplex(0.0,-0.5)*EvtGenFunctions::directProd(dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v),
280 0 : (g1*p4_pi+g2*mb*v))+
281 0 : EvtComplex(0.0,-0.5)*EvtGenFunctions::directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
282 0 : dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v));
283 :
284 0 : l1=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(0));
285 0 : l2=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(1));
286 0 : }
287 : else{
288 0 : report(Severity::Debug,"EvtGen") << "42387dfs878w wrong lepton number\n";
289 : }
290 : }
291 :
292 0 : et0=omega.cont2( dstar->epsParent(0).conj() );
293 0 : et1=omega.cont2( dstar->epsParent(1).conj() );
294 0 : et2=omega.cont2( dstar->epsParent(2).conj() );
295 :
296 0 : vertex(0,0,l1.cont(et0));
297 0 : vertex(0,1,l2.cont(et0));
298 :
299 0 : vertex(1,0,l1.cont(et1));
300 0 : vertex(1,1,l2.cont(et1));
301 :
302 0 : vertex(2,0,l1.cont(et2));
303 0 : vertex(2,1,l2.cont(et2));
304 :
305 : return;
306 :
307 0 : }
308 :
309 : void EvtGoityRoberts::DecayBDpilnuGR(EvtParticle *pb,EvtId nd,
310 : EvtId nlep, EvtId /*nnu*/)
311 :
312 : {
313 : //added by Lange Jan4,2000
314 0 : static EvtId EM=EvtPDL::getId("e-");
315 0 : static EvtId EP=EvtPDL::getId("e+");
316 0 : static EvtId MUM=EvtPDL::getId("mu-");
317 0 : static EvtId MUP=EvtPDL::getId("mu+");
318 :
319 : EvtParticle *d, *pion, *lepton, *neutrino;
320 :
321 0 : pb->initializePhaseSpace(getNDaug(),getDaugs());
322 0 : d=pb->getDaug(0);
323 0 : pion=pb->getDaug(1);
324 0 : lepton=pb->getDaug(2);
325 0 : neutrino=pb->getDaug(3);
326 :
327 0 : EvtVector4C l1, l2, et0, et1, et2;
328 :
329 0 : EvtVector4R v,vp,p4_pi;
330 : double w;
331 :
332 0 : v.set(1.0,0.0,0.0,0.0); //4-velocity of B meson
333 0 : vp=(1.0/d->getP4().mass())*d->getP4(); //4-velocity of D
334 0 : p4_pi=pion->getP4(); //4-momentum of pion
335 0 : w=v*vp; //four velocity transfer.
336 :
337 0 : double mb=EvtPDL::getMeanMass(pb->getId()); //B mass
338 0 : double md=EvtPDL::getMeanMass(nd); //D* mass
339 0 : 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 0 : EvtComplex dmt3(0.563,-0.5*0.191); // Mass splitting = dmt - iGamma/2
351 0 : EvtComplex dmt1(0.392,-0.5*1.040);
352 0 : 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 0 : double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
364 0 : 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 0 : double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
367 0 : pow((2*betas*betap/(betasp)),2.5)*
368 0 : exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
369 0 : double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
370 0 : pow((2*betas*betad/(betasd)),3.5)*
371 0 : exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
372 :
373 0 : EvtComplex h,a1,a2,a3;
374 0 : EvtComplex hnr,a1nr,a2nr,a3nr;
375 0 : EvtComplex hr,a1r,a2r,a3r;
376 :
377 : // Non-resonance part (D* and D** removed by hand - alainb)
378 0 : hnr = g*xi*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb*md);
379 0 : a1nr= -1.0*g*xi*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0);
380 0 : a2nr= g*xi*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb);
381 0 : a3nr=EvtComplex(0.0,0.0);
382 :
383 : // Resonance part (D** remove by hand - alainb)
384 0 : hr = alpha2*rho2*(w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0*mb*md) +
385 0 : alpha3*xi1*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb*md);
386 0 : a1r= -1.0*alpha2*rho2*(w*w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0) -
387 0 : alpha3*xi1*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0);
388 0 : a2r= alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*mb) +
389 0 : alpha2*rho2*(0.5*p4_pi*(w*vp-v)+p4_pi*(vp-w*v))/
390 0 : (3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
391 0 : alpha3*xi1*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb);
392 0 : a3r= -1.0*alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*md) -
393 0 : alpha2*rho2*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmt2))/(2*f0*md);
394 :
395 : // Sum
396 0 : h=hnr+hr;
397 0 : a1=a1nr+a1r;
398 0 : a2=a2nr+a2r;
399 0 : a3=a3nr+a3r;
400 :
401 0 : EvtVector4C omega;
402 :
403 0 : if ( nlep==EM|| nlep==MUM ) {
404 0 : omega=EvtComplex(0.0,-1.0)*h*mb*md*dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v)+
405 0 : a1*p4_pi+a2*mb*v+a3*md*vp;
406 0 : l1=EvtLeptonVACurrent(
407 0 : lepton->spParent(0),neutrino->spParentNeutrino());
408 0 : l2=EvtLeptonVACurrent(
409 0 : lepton->spParent(1),neutrino->spParentNeutrino());
410 0 : }
411 : else{
412 0 : if ( nlep==EP|| nlep==MUP ) {
413 0 : omega=EvtComplex(0.0,1.0)*h*mb*md*dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v)+
414 0 : a1*p4_pi+a2*mb*v+a3*md*vp;
415 0 : l1=EvtLeptonVACurrent(
416 0 : neutrino->spParentNeutrino(),lepton->spParent(0));
417 0 : l2=EvtLeptonVACurrent(
418 0 : neutrino->spParentNeutrino(),lepton->spParent(1));
419 0 : }
420 : else{
421 0 : report(Severity::Error,"EvtGen") << "42387dfs878w wrong lepton number\n";
422 : }
423 : }
424 :
425 0 : vertex(0,l1*omega);
426 0 : vertex(1,l2*omega);
427 :
428 : return;
429 :
430 0 : }
431 :
|
