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: EvtSemiLeptonicBaryonAmp.cc
12 : //
13 : // Description: Routine to implement semileptonic decays of Dirac baryons
14 : //
15 : // Modification history:
16 : //
17 : // R.J. Tesarek May 28, 2004 Module created
18 : // Karen Gibson 1/20/2006 Module updated for 1/2+->1/2+,
19 : // 1/2+->1/2-, 1/2+->3/2- Lambda decays
20 : //
21 : //--------------------------------------------------------------------------
22 :
23 : #include "EvtGenBase/EvtPatches.hh"
24 : #include "EvtGenBase/EvtParticle.hh"
25 : #include "EvtGenBase/EvtGenKine.hh"
26 : #include "EvtGenBase/EvtPDL.hh"
27 : #include "EvtGenBase/EvtReport.hh"
28 : #include "EvtGenBase/EvtTensor4C.hh"
29 : #include "EvtGenBase/EvtVector4C.hh"
30 : #include "EvtGenBase/EvtDiracSpinor.hh"
31 : #include "EvtGenBase/EvtDiracParticle.hh"
32 : #include "EvtGenBase/EvtRaritaSchwinger.hh"
33 : #include "EvtGenBase/EvtSemiLeptonicBaryonAmp.hh"
34 : #include "EvtGenBase/EvtId.hh"
35 : #include "EvtGenBase/EvtAmp.hh"
36 : #include "EvtGenBase/EvtSemiLeptonicFF.hh"
37 : #include "EvtGenBase/EvtGammaMatrix.hh"
38 :
39 : #include <stdlib.h>
40 :
41 : using std::endl;
42 :
43 :
44 0 : EvtSemiLeptonicBaryonAmp::~EvtSemiLeptonicBaryonAmp(){
45 0 : }
46 :
47 :
48 : void EvtSemiLeptonicBaryonAmp::CalcAmp( EvtParticle *parent,
49 : EvtAmp& amp,
50 : EvtSemiLeptonicFF *FormFactors ) {
51 :
52 0 : static EvtId EM=EvtPDL::getId("e-");
53 0 : static EvtId MUM=EvtPDL::getId("mu-");
54 0 : static EvtId TAUM=EvtPDL::getId("tau-");
55 0 : static EvtId EP=EvtPDL::getId("e+");
56 0 : static EvtId MUP=EvtPDL::getId("mu+");
57 0 : static EvtId TAUP=EvtPDL::getId("tau+");
58 :
59 :
60 : //Add the lepton and neutrino 4 momenta to find q2
61 :
62 0 : EvtVector4R q = parent->getDaug(1)->getP4()
63 0 : + parent->getDaug(2)->getP4();
64 0 : double q2 = (q.mass2());
65 :
66 0 : double f1v,f1a,f2v,f2a;
67 0 : double m_meson = parent->getDaug(0)->mass();
68 :
69 0 : FormFactors->getbaryonff(parent->getId(),
70 0 : parent->getDaug(0)->getId(),
71 : q2,
72 : m_meson,
73 : &f1v,
74 : &f1a,
75 : &f2v,
76 : &f2a);
77 :
78 0 : EvtVector4R p4b;
79 0 : p4b.set(parent->mass(),0.0,0.0,0.0);
80 :
81 0 : EvtVector4C temp_00_term1;
82 0 : EvtVector4C temp_00_term2;
83 :
84 0 : EvtVector4C temp_01_term1;
85 0 : EvtVector4C temp_01_term2;
86 :
87 0 : EvtVector4C temp_10_term1;
88 0 : EvtVector4C temp_10_term2;
89 :
90 0 : EvtVector4C temp_11_term1;
91 0 : EvtVector4C temp_11_term2;
92 :
93 0 : EvtDiracSpinor p0=parent->sp(0);
94 0 : EvtDiracSpinor p1=parent->sp(1);
95 :
96 0 : EvtDiracSpinor d0=parent->getDaug(0)->spParent(0);
97 0 : EvtDiracSpinor d1=parent->getDaug(0)->spParent(1);
98 :
99 0 : temp_00_term1.set(0,f1v*(d0*(EvtGammaMatrix::g0()*p0)));
100 0 : temp_00_term2.set(0,f1a*(d0*((EvtGammaMatrix::g0()*EvtGammaMatrix::g5())*p0)));
101 0 : temp_01_term1.set(0,f1v*(d0*(EvtGammaMatrix::g0()*p1)));
102 0 : temp_01_term2.set(0,f1a*(d0*((EvtGammaMatrix::g0()*EvtGammaMatrix::g5())*p1)));
103 0 : temp_10_term1.set(0,f1v*(d1*(EvtGammaMatrix::g0()*p0)));
104 0 : temp_10_term2.set(0,f1a*(d1*((EvtGammaMatrix::g0()*EvtGammaMatrix::g5())*p0)));
105 0 : temp_11_term1.set(0,f1v*(d1*(EvtGammaMatrix::g0()*p1)));
106 0 : temp_11_term2.set(0,f1a*(d1*((EvtGammaMatrix::g0()*EvtGammaMatrix::g5())*p1)));
107 :
108 0 : temp_00_term1.set(1,f1v*(d0*(EvtGammaMatrix::g1()*p0)));
109 0 : temp_00_term2.set(1,f1a*(d0*((EvtGammaMatrix::g1()*EvtGammaMatrix::g5())*p0)));
110 0 : temp_01_term1.set(1,f1v*(d0*(EvtGammaMatrix::g1()*p1)));
111 0 : temp_01_term2.set(1,f1a*(d0*((EvtGammaMatrix::g1()*EvtGammaMatrix::g5())*p1)));
112 0 : temp_10_term1.set(1,f1v*(d1*(EvtGammaMatrix::g1()*p0)));
113 0 : temp_10_term2.set(1,f1a*(d1*((EvtGammaMatrix::g1()*EvtGammaMatrix::g5())*p0)));
114 0 : temp_11_term1.set(1,f1v*(d1*(EvtGammaMatrix::g1()*p1)));
115 0 : temp_11_term2.set(1,f1a*(d1*((EvtGammaMatrix::g1()*EvtGammaMatrix::g5())*p1)));
116 :
117 0 : temp_00_term1.set(2,f1v*(d0*(EvtGammaMatrix::g2()*p0)));
118 0 : temp_00_term2.set(2,f1a*(d0*((EvtGammaMatrix::g2()*EvtGammaMatrix::g5())*p0)));
119 0 : temp_01_term1.set(2,f1v*(d0*(EvtGammaMatrix::g2()*p1)));
120 0 : temp_01_term2.set(2,f1a*(d0*((EvtGammaMatrix::g2()*EvtGammaMatrix::g5())*p1)));
121 0 : temp_10_term1.set(2,f1v*(d1*(EvtGammaMatrix::g2()*p0)));
122 0 : temp_10_term2.set(2,f1a*(d1*((EvtGammaMatrix::g2()*EvtGammaMatrix::g5())*p0)));
123 0 : temp_11_term1.set(2,f1v*(d1*(EvtGammaMatrix::g2()*p1)));
124 0 : temp_11_term2.set(2,f1a*(d1*((EvtGammaMatrix::g2()*EvtGammaMatrix::g5())*p1)));
125 :
126 0 : temp_00_term1.set(3,f1v*(d0*(EvtGammaMatrix::g3()*p0)));
127 0 : temp_00_term2.set(3,f1a*(d0*((EvtGammaMatrix::g3()*EvtGammaMatrix::g5())*p0)));
128 0 : temp_01_term1.set(3,f1v*(d0*(EvtGammaMatrix::g3()*p1)));
129 0 : temp_01_term2.set(3,f1a*(d0*((EvtGammaMatrix::g3()*EvtGammaMatrix::g5())*p1)));
130 0 : temp_10_term1.set(3,f1v*(d1*(EvtGammaMatrix::g3()*p0)));
131 0 : temp_10_term2.set(3,f1a*(d1*((EvtGammaMatrix::g3()*EvtGammaMatrix::g5())*p0)));
132 0 : temp_11_term1.set(3,f1v*(d1*(EvtGammaMatrix::g3()*p1)));
133 0 : temp_11_term2.set(3,f1a*(d1*((EvtGammaMatrix::g3()*EvtGammaMatrix::g5())*p1)));
134 :
135 :
136 :
137 0 : EvtVector4C l1,l2;
138 :
139 0 : EvtId l_num = parent->getDaug(1)->getId();
140 0 : if (l_num==EM||l_num==MUM||l_num==TAUM){
141 :
142 0 : l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
143 0 : parent->getDaug(2)->spParentNeutrino());
144 0 : l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
145 0 : parent->getDaug(2)->spParentNeutrino());
146 0 : }
147 : else{
148 0 : if (l_num==EP||l_num==MUP||l_num==TAUP){
149 0 : l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
150 0 : parent->getDaug(1)->spParent(0));
151 0 : l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
152 0 : parent->getDaug(1)->spParent(1));
153 0 : }
154 : else{
155 0 : report(Severity::Error,"EvtGen") << "Wrong lepton number"<<endl;
156 : }
157 : }
158 :
159 0 : amp.vertex(0,0,0,l1.cont(temp_00_term1+temp_00_term2));
160 0 : amp.vertex(0,0,1,l2.cont(temp_00_term1+temp_00_term2));
161 :
162 0 : amp.vertex(0,1,0,l1.cont(temp_01_term1+temp_01_term2));
163 0 : amp.vertex(0,1,1,l2.cont(temp_01_term1+temp_01_term2));
164 :
165 0 : amp.vertex(1,0,0,l1.cont(temp_10_term1+temp_10_term2));
166 0 : amp.vertex(1,0,1,l2.cont(temp_10_term1+temp_10_term2));
167 :
168 0 : amp.vertex(1,1,0,l1.cont(temp_11_term1+temp_11_term2));
169 0 : amp.vertex(1,1,1,l2.cont(temp_11_term1+temp_11_term2));
170 :
171 : return;
172 0 : }
173 :
174 :
175 :
176 : double EvtSemiLeptonicBaryonAmp::CalcMaxProb( EvtId parent, EvtId baryon,
177 : EvtId lepton, EvtId nudaug,
178 : EvtSemiLeptonicFF *FormFactors,
179 : EvtComplex r00, EvtComplex r01,
180 : EvtComplex r10, EvtComplex r11) {
181 :
182 : //This routine takes the arguements parent, baryon, and lepton
183 : //number, and a form factor model, and returns a maximum
184 : //probability for this semileptonic form factor model. A
185 : //brute force method is used. The 2D cos theta lepton and
186 : //q2 phase space is probed.
187 :
188 : //Start by declaring a particle at rest.
189 :
190 : //It only makes sense to have a scalar parent. For now.
191 : //This should be generalized later.
192 :
193 : // EvtScalarParticle *scalar_part;
194 : // scalar_part=new EvtScalarParticle;
195 :
196 : EvtDiracParticle *dirac_part;
197 : EvtParticle *root_part;
198 0 : dirac_part=new EvtDiracParticle;
199 :
200 : //cludge to avoid generating random numbers!
201 : // scalar_part->noLifeTime();
202 0 : dirac_part->noLifeTime();
203 :
204 0 : EvtVector4R p_init;
205 :
206 0 : p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
207 : // scalar_part->init(parent,p_init);
208 : // root_part=(EvtParticle *)scalar_part;
209 : // root_part->set_type(EvtSpinType::SCALAR);
210 :
211 0 : dirac_part->init(parent,p_init);
212 : root_part=(EvtParticle *)dirac_part;
213 0 : root_part->setDiagonalSpinDensity();
214 :
215 : EvtParticle *daughter, *lep, *trino;
216 :
217 0 : EvtAmp amp;
218 :
219 0 : EvtId listdaug[3];
220 0 : listdaug[0] = baryon;
221 0 : listdaug[1] = lepton;
222 0 : listdaug[2] = nudaug;
223 :
224 0 : amp.init(parent,3,listdaug);
225 :
226 0 : root_part->makeDaughters(3,listdaug);
227 0 : daughter=root_part->getDaug(0);
228 0 : lep=root_part->getDaug(1);
229 0 : trino=root_part->getDaug(2);
230 :
231 : //cludge to avoid generating random numbers!
232 0 : daughter->noLifeTime();
233 0 : lep->noLifeTime();
234 0 : trino->noLifeTime();
235 :
236 :
237 : //Initial particle is unpolarized, well it is a scalar so it is
238 : //trivial
239 0 : EvtSpinDensity rho;
240 0 : rho.setDiag(root_part->getSpinStates());
241 :
242 : double mass[3];
243 :
244 0 : double m = root_part->mass();
245 :
246 0 : EvtVector4R p4baryon, p4lepton, p4nu, p4w;
247 : double q2max;
248 :
249 : double q2, elepton, plepton;
250 : int i,j;
251 : double erho,prho,costl;
252 :
253 : double maxfoundprob = 0.0;
254 : double prob = -10.0;
255 : int massiter;
256 :
257 0 : for (massiter=0;massiter<3;massiter++){
258 :
259 0 : mass[0] = EvtPDL::getMass(baryon);
260 0 : mass[1] = EvtPDL::getMass(lepton);
261 0 : mass[2] = EvtPDL::getMass(nudaug);
262 0 : if ( massiter==1 ) {
263 0 : mass[0] = EvtPDL::getMinMass(baryon);
264 0 : }
265 0 : if ( massiter==2 ) {
266 0 : mass[0] = EvtPDL::getMaxMass(baryon);
267 0 : }
268 :
269 0 : q2max = (m-mass[0])*(m-mass[0]);
270 :
271 : //loop over q2
272 :
273 0 : for (i=0;i<25;i++) {
274 0 : q2 = ((i+0.5)*q2max)/25.0;
275 :
276 0 : erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
277 :
278 0 : prho = sqrt(erho*erho-mass[0]*mass[0]);
279 :
280 0 : p4baryon.set(erho,0.0,0.0,-1.0*prho);
281 0 : p4w.set(m-erho,0.0,0.0,prho);
282 :
283 : //This is in the W rest frame
284 0 : elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
285 0 : plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
286 :
287 0 : double probctl[3];
288 :
289 0 : for (j=0;j<3;j++) {
290 :
291 0 : costl = 0.99*(j - 1.0);
292 :
293 : //These are in the W rest frame. Need to boost out into
294 : //the B frame.
295 0 : p4lepton.set(elepton,0.0,
296 0 : plepton*sqrt(1.0-costl*costl),plepton*costl);
297 0 : p4nu.set(plepton,0.0,
298 0 : -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
299 :
300 0 : EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
301 0 : p4lepton=boostTo(p4lepton,boost);
302 0 : p4nu=boostTo(p4nu,boost);
303 :
304 : //Now initialize the daughters...
305 :
306 0 : daughter->init(baryon,p4baryon);
307 0 : lep->init(lepton,p4lepton);
308 0 : trino->init(nudaug,p4nu);
309 :
310 0 : CalcAmp(root_part,amp,FormFactors,r00,r01,r10,r11);
311 :
312 : //Now find the probability at this q2 and cos theta lepton point
313 : //and compare to maxfoundprob.
314 :
315 : //Do a little magic to get the probability!!
316 0 : prob = rho.normalizedProb(amp.getSpinDensity());
317 :
318 0 : probctl[j]=prob;
319 0 : }
320 :
321 : //probclt contains prob at ctl=-1,0,1.
322 : //prob=a+b*ctl+c*ctl^2
323 :
324 0 : double a=probctl[1];
325 0 : double b=0.5*(probctl[2]-probctl[0]);
326 0 : double c=0.5*(probctl[2]+probctl[0])-probctl[1];
327 :
328 : prob=probctl[0];
329 0 : if (probctl[1]>prob) prob=probctl[1];
330 0 : if (probctl[2]>prob) prob=probctl[2];
331 :
332 0 : if (fabs(c)>1e-20){
333 0 : double ctlx=-0.5*b/c;
334 0 : if (fabs(ctlx)<1.0){
335 0 : double probtmp=a+b*ctlx+c*ctlx*ctlx;
336 0 : if (probtmp>prob) prob=probtmp;
337 0 : }
338 :
339 0 : }
340 :
341 : //report(Severity::Debug,"EvtGen") << "prob,probctl:"<<prob<<" "
342 : // << probctl[0]<<" "
343 : // << probctl[1]<<" "
344 : // << probctl[2]<<std::endl;
345 :
346 0 : if ( prob > maxfoundprob ) {
347 : maxfoundprob = prob;
348 0 : }
349 :
350 0 : }
351 0 : if ( EvtPDL::getWidth(baryon) <= 0.0 ) {
352 : //if the particle is narrow dont bother with changing the mass.
353 : massiter = 4;
354 0 : }
355 :
356 : }
357 0 : root_part->deleteTree();
358 :
359 0 : maxfoundprob *=1.1;
360 : return maxfoundprob;
361 :
362 0 : }
363 : void EvtSemiLeptonicBaryonAmp::CalcAmp(EvtParticle *parent,
364 : EvtAmp& amp,
365 : EvtSemiLeptonicFF *FormFactors,
366 : EvtComplex r00, EvtComplex r01,
367 : EvtComplex r10, EvtComplex r11) {
368 : // Leptons
369 0 : static EvtId EM=EvtPDL::getId("e-");
370 0 : static EvtId MUM=EvtPDL::getId("mu-");
371 0 : static EvtId TAUM=EvtPDL::getId("tau-");
372 : // Anti-Leptons
373 0 : static EvtId EP=EvtPDL::getId("e+");
374 0 : static EvtId MUP=EvtPDL::getId("mu+");
375 0 : static EvtId TAUP=EvtPDL::getId("tau+");
376 :
377 : // Baryons
378 0 : static EvtId LAMCP=EvtPDL::getId("Lambda_c+");
379 0 : static EvtId LAMC1P=EvtPDL::getId("Lambda_c(2593)+");
380 0 : static EvtId LAMC2P=EvtPDL::getId("Lambda_c(2625)+");
381 0 : static EvtId LAMB=EvtPDL::getId("Lambda_b0");
382 :
383 : // Anti-Baryons
384 0 : static EvtId LAMCM=EvtPDL::getId("anti-Lambda_c-");
385 0 : static EvtId LAMC1M=EvtPDL::getId("anti-Lambda_c(2593)-");
386 0 : static EvtId LAMC2M=EvtPDL::getId("anti-Lambda_c(2625)-");
387 0 : static EvtId LAMBB=EvtPDL::getId("anti-Lambda_b0");
388 :
389 : // Set the spin density matrix of the parent baryon
390 0 : EvtSpinDensity rho;
391 0 : rho.setDim(2);
392 0 : rho.set(0,0,r00);
393 0 : rho.set(0,1,r01);
394 :
395 0 : rho.set(1,0,r10);
396 0 : rho.set(1,1,r11);
397 :
398 0 : EvtVector4R vector4P = parent->getP4Lab();
399 0 : double pmag = vector4P.d3mag();
400 0 : double cosTheta = vector4P.get(3)/pmag;
401 :
402 0 : double theta = acos(cosTheta);
403 0 : double phi = atan2(vector4P.get(2), vector4P.get(1));
404 :
405 0 : parent->setSpinDensityForwardHelicityBasis(rho,phi,theta, 0.0);
406 : //parent->setSpinDensityForward(rho);
407 :
408 : // Set the four momentum of the parent baryon in it's rest frame
409 0 : EvtVector4R p4b;
410 0 : p4b.set(parent->mass(), 0.0,0.0,0.0);
411 :
412 : // Get the four momentum of the daughter baryon in the parent's rest frame
413 0 : EvtVector4R p4daught = parent->getDaug(0)->getP4();
414 :
415 : // Add the lepton and neutrino 4 momenta to find q (q^2)
416 0 : EvtVector4R q = parent->getDaug(1)->getP4()
417 0 : + parent->getDaug(2)->getP4();
418 :
419 0 : double q2 = q.mass2();
420 :
421 :
422 0 : EvtId l_num = parent->getDaug(1)->getId();
423 0 : EvtId bar_num = parent->getDaug(0)->getId();
424 0 : EvtId par_num = parent->getId();
425 :
426 0 : double baryonmass = parent->getDaug(0)->mass();
427 :
428 : // Handle spin-1/2 daughter baryon Dirac spinor cases
429 0 : if( EvtPDL::getSpinType(parent->getDaug(0)->getId())==EvtSpinType::DIRAC ) {
430 :
431 : // Set the form factors
432 0 : double f1,f2,f3,g1,g2,g3;
433 0 : FormFactors->getdiracff( par_num,
434 0 : bar_num,
435 : q2,
436 : baryonmass,
437 : &f1, &f2, &f3,
438 : &g1, &g2, &g3);
439 :
440 0 : const double form_fact[6] = {f1, f2, f3, g1, g2, g3};
441 :
442 0 : EvtVector4C b11, b12, b21, b22, l1, l2;
443 :
444 : // Lepton Current
445 0 : if(l_num==EM || l_num==MUM || l_num==TAUM){
446 :
447 0 : l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
448 0 : parent->getDaug(2)->spParentNeutrino());
449 0 : l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
450 0 : parent->getDaug(2)->spParentNeutrino());
451 :
452 0 : } else if (l_num==EP || l_num==MUP || l_num==TAUP) {
453 :
454 0 : l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
455 0 : parent->getDaug(1)->spParent(0));
456 0 : l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
457 0 : parent->getDaug(1)->spParent(1));
458 :
459 : } else {
460 0 : report(Severity::Error,"EvtGen")<< "Wrong lepton number \n";
461 0 : ::abort();
462 : }
463 :
464 : // Baryon current
465 :
466 : // Flag for particle/anti-particle parent, daughter with same/opp. parity
467 : // pflag = 0 => particle, same parity parent, daughter
468 : // pflag = 1 => particle, opp. parity parent, daughter
469 : // pflag = 2 => anti-particle, same parity parent, daughter
470 : // pflag = 3 => anti-particle, opp. parity parent, daughter
471 :
472 : int pflag = 0;
473 :
474 : // Handle 1/2+ -> 1/2+ first
475 0 : if ( (par_num==LAMB && bar_num==LAMCP)
476 0 : || (par_num==LAMBB && bar_num==LAMCM) ) {
477 :
478 : // Set particle/anti-particle flag
479 0 : if (bar_num==LAMCP)
480 0 : pflag = 0;
481 0 : else if (bar_num==LAMCM)
482 0 : pflag = 2;
483 :
484 0 : b11=EvtBaryonVACurrent(parent->getDaug(0)->spParent(0),
485 0 : parent->sp(0),
486 0 : p4b, p4daught, form_fact, pflag);
487 0 : b21=EvtBaryonVACurrent(parent->getDaug(0)->spParent(0),
488 0 : parent->sp(1),
489 0 : p4b, p4daught, form_fact, pflag);
490 0 : b12=EvtBaryonVACurrent(parent->getDaug(0)->spParent(1),
491 0 : parent->sp(0),
492 0 : p4b, p4daught, form_fact, pflag);
493 0 : b22=EvtBaryonVACurrent(parent->getDaug(0)->spParent(1),
494 0 : parent->sp(1),
495 0 : p4b, p4daught, form_fact, pflag);
496 0 : }
497 :
498 : // Handle 1/2+ -> 1/2- second
499 0 : else if( (par_num==LAMB && bar_num==LAMC1P)
500 0 : || (par_num==LAMBB && bar_num==LAMC1M) ) {
501 :
502 : // Set particle/anti-particle flag
503 0 : if (bar_num==LAMC1P)
504 0 : pflag = 1;
505 0 : else if (bar_num==LAMC1M)
506 0 : pflag = 3;
507 :
508 0 : b11=EvtBaryonVACurrent((parent->getDaug(0)->spParent(0)),
509 0 : (EvtGammaMatrix::g5()*parent->sp(0)),
510 0 : p4b, p4daught, form_fact, pflag);
511 0 : b21=EvtBaryonVACurrent((parent->getDaug(0)->spParent(0)),
512 0 : (EvtGammaMatrix::g5()*parent->sp(1)),
513 0 : p4b, p4daught, form_fact, pflag);
514 0 : b12=EvtBaryonVACurrent((parent->getDaug(0)->spParent(1)),
515 0 : (EvtGammaMatrix::g5()*parent->sp(0)),
516 0 : p4b, p4daught, form_fact, pflag);
517 0 : b22=EvtBaryonVACurrent((parent->getDaug(0)->spParent(1)),
518 0 : (EvtGammaMatrix::g5()*parent->sp(1)),
519 0 : p4b, p4daught, form_fact, pflag);
520 :
521 : }
522 :
523 : else {
524 0 : report(Severity::Error,"EvtGen") << "Dirac semilep. baryon current "
525 0 : << "not implemented for this decay sequence."
526 0 : << std::endl;
527 0 : ::abort();
528 : }
529 :
530 0 : amp.vertex(0,0,0,l1*b11);
531 0 : amp.vertex(0,0,1,l2*b11);
532 :
533 0 : amp.vertex(1,0,0,l1*b21);
534 0 : amp.vertex(1,0,1,l2*b21);
535 :
536 0 : amp.vertex(0,1,0,l1*b12);
537 0 : amp.vertex(0,1,1,l2*b12);
538 :
539 0 : amp.vertex(1,1,0,l1*b22);
540 0 : amp.vertex(1,1,1,l2*b22);
541 :
542 0 : }
543 :
544 : // Need special handling for the spin-3/2 daughter baryon
545 : // Rarita-Schwinger spinor cases
546 0 : else if( EvtPDL::getSpinType(parent->getDaug(0)->getId())==EvtSpinType::RARITASCHWINGER ) {
547 :
548 : // Set the form factors
549 0 : double f1,f2,f3,f4,g1,g2,g3,g4;
550 0 : FormFactors->getraritaff( par_num,
551 0 : bar_num,
552 : q2,
553 : baryonmass,
554 : &f1, &f2, &f3, &f4,
555 : &g1, &g2, &g3, &g4);
556 :
557 0 : const double form_fact[8] = {f1, f2, f3, f4, g1, g2, g3, g4};
558 :
559 0 : EvtId l_num = parent->getDaug(1)->getId();
560 :
561 0 : EvtVector4C b11, b12, b21, b22, b13, b23, b14, b24, l1, l2;
562 :
563 : // Lepton Current
564 0 : if (l_num==EM || l_num==MUM || l_num==TAUM) {
565 : // Lepton Current
566 0 : l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
567 0 : parent->getDaug(2)->spParentNeutrino());
568 0 : l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
569 0 : parent->getDaug(2)->spParentNeutrino());
570 0 : }
571 0 : else if (l_num==EP || l_num==MUP || l_num==TAUP) {
572 0 : l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
573 0 : parent->getDaug(1)->spParent(0));
574 0 : l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
575 0 : parent->getDaug(1)->spParent(1));
576 :
577 0 : } else {
578 0 : report(Severity::Error,"EvtGen")<< "Wrong lepton number \n";
579 : }
580 :
581 : // Baryon Current
582 : // Declare particle, anti-particle flag, same/opp. parity
583 : // pflag = 0 => particle
584 : // pflag = 1 => anti-particle
585 : int pflag = 0;
586 :
587 : // Handle cases of 1/2+ -> 3/2-
588 0 : if (par_num==LAMB && bar_num==LAMC2P) {
589 : // Set flag for particle case
590 : pflag = 0;
591 0 : }
592 0 : else if (par_num==LAMBB && bar_num==LAMC2M) {
593 : // Set flag for anti-particle case
594 : pflag = 1;
595 : }
596 : else {
597 0 : report(Severity::Error,"EvtGen") << "Rarita-Schwinger semilep. baryon current "
598 0 : << "not implemented for this decay sequence."
599 0 : << std::endl;
600 0 : ::abort();
601 : }
602 :
603 : // Baryon current
604 0 : b11=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(0),
605 0 : parent->sp(0),
606 0 : p4b, p4daught, form_fact, pflag);
607 0 : b21=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(0),
608 0 : parent->sp(1),
609 0 : p4b, p4daught, form_fact, pflag);
610 :
611 0 : b12=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(1),
612 0 : parent->sp(0),
613 0 : p4b, p4daught, form_fact, pflag);
614 0 : b22=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(1),
615 0 : parent->sp(1),
616 0 : p4b, p4daught, form_fact, pflag);
617 :
618 0 : b13=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(2),
619 0 : parent->sp(0),
620 0 : p4b, p4daught, form_fact, pflag);
621 0 : b23=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(2),
622 0 : parent->sp(1),
623 0 : p4b, p4daught, form_fact, pflag);
624 :
625 0 : b14=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(3),
626 0 : parent->sp(0),
627 0 : p4b, p4daught, form_fact, pflag);
628 0 : b24=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(3),
629 0 : parent->sp(1),
630 0 : p4b, p4daught, form_fact, pflag);
631 :
632 0 : amp.vertex(0,0,0,l1*b11);
633 0 : amp.vertex(0,0,1,l2*b11);
634 :
635 0 : amp.vertex(1,0,0,l1*b21);
636 0 : amp.vertex(1,0,1,l2*b21);
637 :
638 0 : amp.vertex(0,1,0,l1*b12);
639 0 : amp.vertex(0,1,1,l2*b12);
640 :
641 0 : amp.vertex(1,1,0,l1*b22);
642 0 : amp.vertex(1,1,1,l2*b22);
643 :
644 0 : amp.vertex(0,2,0,l1*b13);
645 0 : amp.vertex(0,2,1,l2*b13);
646 :
647 0 : amp.vertex(1,2,0,l1*b23);
648 0 : amp.vertex(1,2,1,l2*b23);
649 :
650 0 : amp.vertex(0,3,0,l1*b14);
651 0 : amp.vertex(0,3,1,l2*b14);
652 :
653 0 : amp.vertex(1,3,0,l1*b24);
654 0 : amp.vertex(1,3,1,l2*b24);
655 :
656 0 : }
657 :
658 0 : }
659 :
660 :
661 : EvtVector4C EvtSemiLeptonicBaryonAmp::EvtBaryonVACurrent( const EvtDiracSpinor& Bf,
662 : const EvtDiracSpinor& Bi,
663 : EvtVector4R parent,
664 : EvtVector4R daught,
665 : const double *ff,
666 : int pflag) {
667 :
668 : // flag == 0 => particle, same parity
669 : // flag == 1 => particle, opposite parity
670 : // flag == 2 => anti-particle, same parity
671 : // flag == 3 => anti-particle, opposite parity
672 :
673 : // particle
674 0 : EvtComplex cv = EvtComplex(1.0, 0.);
675 0 : EvtComplex ca = EvtComplex(1.0, 0.);
676 :
677 0 : EvtComplex cg0 = EvtComplex(1.0, 0.);
678 0 : EvtComplex cg5 = EvtComplex(1.0, 0.);
679 :
680 : // antiparticle- same parity parent & daughter
681 0 : if( pflag == 2 ) {
682 0 : cv = EvtComplex(-1.0, 0.);
683 0 : ca = EvtComplex(1.0, 0.);
684 :
685 0 : cg0 = EvtComplex(1.0, 0.0);
686 0 : cg5 = EvtComplex(0.0, -1.0);
687 0 : }
688 : // antiparticle- opposite parity parent & daughter
689 0 : else if( pflag == 3) {
690 0 : cv = EvtComplex(1.0, 0.);
691 0 : ca = EvtComplex(-1.0, 0.);
692 :
693 0 : cg0 = EvtComplex(0.0, -1.0);
694 0 : cg5 = EvtComplex(1.0, 0.0);
695 0 : }
696 :
697 0 : EvtVector4C t[6];
698 :
699 : // Term 1 = \bar{u}(p',s')*(F_1(q^2)*\gamma_{mu})*u(p,s)
700 0 : t[0] = cv*EvtLeptonVCurrent( Bf, Bi);
701 :
702 : // Term 2 = \bar{u}(p',s')*(F_2(q^2)*(p_{mu}/m_{\Lambda_Q}))*u(p,s)
703 0 : t[1] = cg0*EvtLeptonSCurrent( Bf, Bi ) * (parent/parent.mass());
704 :
705 : // Term 3 = \bar{u}(p',s')*(F_3(q^2)*(p'_{mu}/m_{\Lambda_q}))*u(p,s)
706 0 : t[2] = cg0*EvtLeptonSCurrent( Bf, Bi ) * (daught/daught.mass());
707 :
708 : // Term 4 = \bar{u}(p',s')*(G_1(q^2)*\gamma_{mu}*\gamma_5)*u(p,s)
709 0 : t[3] = ca*EvtLeptonACurrent( Bf, Bi);
710 :
711 : // Term 5 = \bar{u}(p',s')*(G_2(q^2)*(p_{mu}/m_{\Lambda_Q})*\gamma_5)*u(p,s)
712 0 : t[4] = cg5*EvtLeptonPCurrent( Bf, Bi ) * (parent/parent.mass());
713 :
714 : // Term 6 = \bar{u}(p',s')*(G_3(q^2)*(p'_{mu}/m_{\Lambda_q})*\gamma_5)*u(p,s)
715 0 : t[5] = cg5*EvtLeptonPCurrent( Bf, Bi ) * (daught/daught.mass());
716 :
717 : // Sum the individual terms
718 0 : EvtVector4C current = (ff[0]*t[0] + ff[1]*t[1] + ff[2]*t[2]
719 0 : - ff[3]*t[3] - ff[4]*t[4] - ff[5]*t[5]);
720 :
721 : return current;
722 0 : }
723 :
724 : EvtVector4C EvtSemiLeptonicBaryonAmp::EvtBaryonVARaritaCurrent( const EvtRaritaSchwinger& Bf,
725 : const EvtDiracSpinor& Bi,
726 : EvtVector4R parent,
727 : EvtVector4R daught,
728 : const double *ff,
729 : int pflag) {
730 :
731 : // flag == 0 => particle
732 : // flag == 1 => anti-particle
733 :
734 : // particle
735 0 : EvtComplex cv = EvtComplex(1.0, 0.);
736 0 : EvtComplex ca = EvtComplex(1.0, 0.);
737 :
738 0 : EvtComplex cg0 = EvtComplex(1.0, 0.);
739 0 : EvtComplex cg5 = EvtComplex(1.0, 0.);
740 :
741 : // antiparticle
742 0 : if( pflag == 1 ) {
743 0 : cv = EvtComplex(-1.0, 0.);
744 0 : ca = EvtComplex(1.0, 0.);
745 :
746 0 : cg0 = EvtComplex(1.0, 0.0);
747 0 : cg5 = EvtComplex(0.0, -1.0);
748 0 : }
749 :
750 0 : EvtVector4C t[8];
751 0 : EvtTensor4C id;
752 0 : id.setdiag(1.0,1.0,1.0,1.0);
753 :
754 0 : EvtDiracSpinor tmp;
755 0 : for(int i=0;i<4;i++){
756 0 : tmp.set_spinor(i,Bf.getVector(i)*parent);
757 : }
758 :
759 0 : EvtVector4C v1,v2;
760 0 : for(int i=0;i<4;i++){
761 0 : v1.set(i,EvtLeptonSCurrent(Bf.getSpinor(i),Bi));
762 0 : v2.set(i,EvtLeptonPCurrent(Bf.getSpinor(i),Bi));
763 : }
764 :
765 : // Term 1 = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_1(q^2)*\gamma_{mu})*u(p,s)
766 0 : t[0] = (cv/parent.mass()) * EvtLeptonVCurrent(tmp, Bi);
767 :
768 : // Term 2
769 : // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_2(q^2)*(p_{mu}/m_{\Lambda_Q}))*u(p,s)
770 0 : t[1] = ((cg0/parent.mass()) * EvtLeptonSCurrent(tmp, Bi)) * (parent/parent.mass());
771 :
772 : // Term 3
773 : // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_3(q^2)*(p'_{mu}/m_{\Lambda_q}))*u(p,s)
774 0 : t[2] = ((cg0/parent.mass()) * EvtLeptonSCurrent(tmp, Bi)) * (daught/daught.mass());
775 :
776 : // Term 4 = \bar{u}^{\alpha}(p',s')*(F_4(q^2)*g_{\alpha,\mu})*u(p,s)
777 0 : t[3] = cg0*(id.cont2(v1));
778 :
779 : // Term 5
780 : // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(G_1(q^2)*\gamma_{mu}*\gamma_5)*u(p,s)
781 0 : t[4] = (ca/parent.mass()) * EvtLeptonACurrent(tmp, Bi);
782 :
783 : // Term 6
784 : // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})
785 : // *(G_2(q^2)*(p_{mu}/m_{\Lambda_Q})*\gamma_5)*u(p,s)
786 0 : t[5] = ((cg5/parent.mass()) * EvtLeptonPCurrent(tmp, Bi)) * (parent/parent.mass());
787 :
788 : // Term 7
789 : // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})
790 : // *(G_3(q^2)*(p'_{mu}/m_{\Lambda_q})*\gamma_5)*u(p,s)
791 0 : t[6] = ((cg5/parent.mass()) * EvtLeptonPCurrent(tmp, Bi)) * (daught/daught.mass());
792 :
793 : // Term 8 = \bar{u}^{\alpha}(p',s')*(G_4(q^2)*g_{\alpha,\mu}*\gamma_5))*u(p,s)
794 0 : t[7] = cg5*(id.cont2(v2));
795 :
796 : // Sum the individual terms
797 0 : EvtVector4C current = (ff[0]*t[0] + ff[1]*t[1] + ff[2]*t[2] + ff[3]*t[3]
798 0 : - ff[4]*t[4] - ff[5]*t[5] - ff[6]*t[6] - ff[7]*t[7]);
799 :
800 : return current;
801 0 : }
|
