TEvtGen/EvtGenBase/EvtSemiLeptonicBaryonAmp.cpp

   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 EvtSemiLeptonicBaryonAmp::~EvtSemiLeptonicBaryonAmp(){
  45 }
  46
  47
  48 void EvtSemiLeptonicBaryonAmp::CalcAmp( EvtParticle *parent,
  49                                         EvtAmp& amp,
  50                                         EvtSemiLeptonicFF *FormFactors ) {
  51
  52   static EvtId EM=EvtPDL::getId("e-");
  53   static EvtId MUM=EvtPDL::getId("mu-");
  54   static EvtId TAUM=EvtPDL::getId("tau-");
  55   static EvtId EP=EvtPDL::getId("e+");
  56   static EvtId MUP=EvtPDL::getId("mu+");
  57   static EvtId TAUP=EvtPDL::getId("tau+");
  58
  59
  60   //Add the lepton and neutrino 4 momenta to find q2
  61
  62   EvtVector4R q = parent->getDaug(1)->getP4()
  63                     + parent->getDaug(2)->getP4();
  64   double q2 = (q.mass2());
  65
  66   double f1v,f1a,f2v,f2a;
  67   double m_meson = parent->getDaug(0)->mass();
  68
  69   FormFactors->getbaryonff(parent->getId(),
  70                            parent->getDaug(0)->getId(),
  71                            q2,
  72                            m_meson,
  73                            &f1v,
  74                            &f1a,
  75                            &f2v,
  76                            &f2a);
  77
  78   EvtVector4R p4b;
  79   p4b.set(parent->mass(),0.0,0.0,0.0);
  80
  81   EvtVector4C temp_00_term1;
  82   EvtVector4C temp_00_term2;
  83
  84   EvtVector4C temp_01_term1;
  85   EvtVector4C temp_01_term2;
  86
  87   EvtVector4C temp_10_term1;
  88   EvtVector4C temp_10_term2;
  89
  90   EvtVector4C temp_11_term1;
  91   EvtVector4C temp_11_term2;
  92
  93   EvtDiracSpinor p0=parent->sp(0);
  94   EvtDiracSpinor p1=parent->sp(1);
  95
  96   EvtDiracSpinor d0=parent->getDaug(0)->spParent(0);
  97   EvtDiracSpinor d1=parent->getDaug(0)->spParent(1);
  98
  99   temp_00_term1.set(0,f1v*(d0*(EvtGammaMatrix::g0()*p0)));
 100   temp_00_term2.set(0,f1a*(d0*((EvtGammaMatrix::g0()*EvtGammaMatrix::g5())*p0)));
 101   temp_01_term1.set(0,f1v*(d0*(EvtGammaMatrix::g0()*p1)));
 102   temp_01_term2.set(0,f1a*(d0*((EvtGammaMatrix::g0()*EvtGammaMatrix::g5())*p1)));
 103   temp_10_term1.set(0,f1v*(d1*(EvtGammaMatrix::g0()*p0)));
 104   temp_10_term2.set(0,f1a*(d1*((EvtGammaMatrix::g0()*EvtGammaMatrix::g5())*p0)));
 105   temp_11_term1.set(0,f1v*(d1*(EvtGammaMatrix::g0()*p1)));
 106   temp_11_term2.set(0,f1a*(d1*((EvtGammaMatrix::g0()*EvtGammaMatrix::g5())*p1)));
 107
 108   temp_00_term1.set(1,f1v*(d0*(EvtGammaMatrix::g1()*p0)));
 109   temp_00_term2.set(1,f1a*(d0*((EvtGammaMatrix::g1()*EvtGammaMatrix::g5())*p0)));
 110   temp_01_term1.set(1,f1v*(d0*(EvtGammaMatrix::g1()*p1)));
 111   temp_01_term2.set(1,f1a*(d0*((EvtGammaMatrix::g1()*EvtGammaMatrix::g5())*p1)));
 112   temp_10_term1.set(1,f1v*(d1*(EvtGammaMatrix::g1()*p0)));
 113   temp_10_term2.set(1,f1a*(d1*((EvtGammaMatrix::g1()*EvtGammaMatrix::g5())*p0)));
 114   temp_11_term1.set(1,f1v*(d1*(EvtGammaMatrix::g1()*p1)));
 115   temp_11_term2.set(1,f1a*(d1*((EvtGammaMatrix::g1()*EvtGammaMatrix::g5())*p1)));
 116
 117   temp_00_term1.set(2,f1v*(d0*(EvtGammaMatrix::g2()*p0)));
 118   temp_00_term2.set(2,f1a*(d0*((EvtGammaMatrix::g2()*EvtGammaMatrix::g5())*p0)));
 119   temp_01_term1.set(2,f1v*(d0*(EvtGammaMatrix::g2()*p1)));
 120   temp_01_term2.set(2,f1a*(d0*((EvtGammaMatrix::g2()*EvtGammaMatrix::g5())*p1)));
 121   temp_10_term1.set(2,f1v*(d1*(EvtGammaMatrix::g2()*p0)));
 122   temp_10_term2.set(2,f1a*(d1*((EvtGammaMatrix::g2()*EvtGammaMatrix::g5())*p0)));
 123   temp_11_term1.set(2,f1v*(d1*(EvtGammaMatrix::g2()*p1)));
 124   temp_11_term2.set(2,f1a*(d1*((EvtGammaMatrix::g2()*EvtGammaMatrix::g5())*p1)));
 125
 126   temp_00_term1.set(3,f1v*(d0*(EvtGammaMatrix::g3()*p0)));
 127   temp_00_term2.set(3,f1a*(d0*((EvtGammaMatrix::g3()*EvtGammaMatrix::g5())*p0)));
 128   temp_01_term1.set(3,f1v*(d0*(EvtGammaMatrix::g3()*p1)));
 129   temp_01_term2.set(3,f1a*(d0*((EvtGammaMatrix::g3()*EvtGammaMatrix::g5())*p1)));
 130   temp_10_term1.set(3,f1v*(d1*(EvtGammaMatrix::g3()*p0)));
 131   temp_10_term2.set(3,f1a*(d1*((EvtGammaMatrix::g3()*EvtGammaMatrix::g5())*p0)));
 132   temp_11_term1.set(3,f1v*(d1*(EvtGammaMatrix::g3()*p1)));
 133   temp_11_term2.set(3,f1a*(d1*((EvtGammaMatrix::g3()*EvtGammaMatrix::g5())*p1)));
 134
 135
 136
 137   EvtVector4C l1,l2;
 138
 139   EvtId l_num = parent->getDaug(1)->getId();
 140   if (l_num==EM||l_num==MUM||l_num==TAUM){
 141
 142     l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
 143                           parent->getDaug(2)->spParentNeutrino());
 144     l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
 145                           parent->getDaug(2)->spParentNeutrino());
 146   }
 147   else{
 148     if (l_num==EP||l_num==MUP||l_num==TAUP){
 149     l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
 150                             parent->getDaug(1)->spParent(0));
 151     l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
 152                             parent->getDaug(1)->spParent(1));
 153     }
 154     else{
 155       report(ERROR,"EvtGen") << "Wrong lepton number"<<endl;
 156     }
 157   }
 158
 159   amp.vertex(0,0,0,l1.cont(temp_00_term1+temp_00_term2));
 160   amp.vertex(0,0,1,l2.cont(temp_00_term1+temp_00_term2));
 161
 162   amp.vertex(0,1,0,l1.cont(temp_01_term1+temp_01_term2));
 163   amp.vertex(0,1,1,l2.cont(temp_01_term1+temp_01_term2));
 164
 165   amp.vertex(1,0,0,l1.cont(temp_10_term1+temp_10_term2));
 166   amp.vertex(1,0,1,l2.cont(temp_10_term1+temp_10_term2));
 167
 168   amp.vertex(1,1,0,l1.cont(temp_11_term1+temp_11_term2));
 169   amp.vertex(1,1,1,l2.cont(temp_11_term1+temp_11_term2));
 170
 171   return;
 172 }
 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   dirac_part=new EvtDiracParticle;
 199
 200   //cludge to avoid generating random numbers!
 201   //  scalar_part->noLifeTime();
 202   dirac_part->noLifeTime();
 203
 204   EvtVector4R p_init;
 205
 206   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   dirac_part->init(parent,p_init);
 212   root_part=(EvtParticle *)dirac_part;
 213   root_part->setDiagonalSpinDensity();
 214
 215   EvtParticle *daughter, *lep, *trino;
 216
 217   EvtAmp amp;
 218
 219   EvtId listdaug[3];
 220   listdaug[0] = baryon;
 221   listdaug[1] = lepton;
 222   listdaug[2] = nudaug;
 223
 224   amp.init(parent,3,listdaug);
 225
 226   root_part->makeDaughters(3,listdaug);
 227   daughter=root_part->getDaug(0);
 228   lep=root_part->getDaug(1);
 229   trino=root_part->getDaug(2);
 230
 231   //cludge to avoid generating random numbers!
 232   daughter->noLifeTime();
 233   lep->noLifeTime();
 234   trino->noLifeTime();
 235
 236
 237   //Initial particle is unpolarized, well it is a scalar so it is
 238   //trivial
 239   EvtSpinDensity rho;
 240   rho.setDiag(root_part->getSpinStates());
 241
 242   double mass[3];
 243
 244   double m = root_part->mass();
 245
 246   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   for (massiter=0;massiter<3;massiter++){
 258
 259     mass[0] = EvtPDL::getMass(baryon);
 260     mass[1] = EvtPDL::getMass(lepton);
 261     mass[2] = EvtPDL::getMass(nudaug);
 262     if ( massiter==1 ) {
 263       mass[0] = EvtPDL::getMinMass(baryon);
 264     }
 265     if ( massiter==2 ) {
 266       mass[0] = EvtPDL::getMaxMass(baryon);
 267     }
 268
 269     q2max = (m-mass[0])*(m-mass[0]);
 270
 271     //loop over q2
 272
 273     for (i=0;i<25;i++) {
 274       q2 = ((i+0.5)*q2max)/25.0;
 275
 276       erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
 277
 278       prho = sqrt(erho*erho-mass[0]*mass[0]);
 279
 280       p4baryon.set(erho,0.0,0.0,-1.0*prho);
 281       p4w.set(m-erho,0.0,0.0,prho);
 282
 283       //This is in the W rest frame
 284       elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
 285       plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
 286
 287       double probctl[3];
 288
 289       for (j=0;j<3;j++) {
 290
 291         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         p4lepton.set(elepton,0.0,
 296                   plepton*sqrt(1.0-costl*costl),plepton*costl);
 297         p4nu.set(plepton,0.0,
 298                  -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
 299
 300         EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
 301         p4lepton=boostTo(p4lepton,boost);
 302         p4nu=boostTo(p4nu,boost);
 303
 304         //Now initialize the daughters...
 305
 306         daughter->init(baryon,p4baryon);
 307         lep->init(lepton,p4lepton);
 308         trino->init(nudaug,p4nu);
 309
 310         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         prob = rho.normalizedProb(amp.getSpinDensity());
 317
 318         probctl[j]=prob;
 319       }
 320
 321       //probclt contains prob at ctl=-1,0,1.
 322       //prob=a+b*ctl+c*ctl^2
 323
 324       double a=probctl[1];
 325       double b=0.5*(probctl[2]-probctl[0]);
 326       double c=0.5*(probctl[2]+probctl[0])-probctl[1];
 327
 328       prob=probctl[0];
 329       if (probctl[1]>prob) prob=probctl[1];
 330       if (probctl[2]>prob) prob=probctl[2];
 331
 332       if (fabs(c)>1e-20){
 333         double ctlx=-0.5*b/c;
 334         if (fabs(ctlx)<1.0){
 335           double probtmp=a+b*ctlx+c*ctlx*ctlx;
 336           if (probtmp>prob) prob=probtmp;
 337         }
 338
 339       }
 340
 341       //report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
 342       //                            << probctl[0]<<" "
 343       //                            << probctl[1]<<" "
 344       //                            << probctl[2]<<std::endl;
 345
 346       if ( prob > maxfoundprob ) {
 347         maxfoundprob = prob;
 348       }
 349
 350     }
 351     if ( EvtPDL::getWidth(baryon) <= 0.0 ) {
 352       //if the particle is narrow dont bother with changing the mass.
 353       massiter = 4;
 354     }
 355
 356   }
 357   root_part->deleteTree();
 358
 359   maxfoundprob *=1.1;
 360   return maxfoundprob;
 361
 362 }
 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   static EvtId EM=EvtPDL::getId("e-");
 370   static EvtId MUM=EvtPDL::getId("mu-");
 371   static EvtId TAUM=EvtPDL::getId("tau-");
 372   //  Anti-Leptons
 373   static EvtId EP=EvtPDL::getId("e+");
 374   static EvtId MUP=EvtPDL::getId("mu+");
 375   static EvtId TAUP=EvtPDL::getId("tau+");
 376
 377   //  Baryons
 378   static EvtId LAMCP=EvtPDL::getId("Lambda_c+");
 379   static EvtId LAMC1P=EvtPDL::getId("Lambda_c(2593)+");
 380   static EvtId LAMC2P=EvtPDL::getId("Lambda_c(2625)+");
 381   static EvtId LAMB=EvtPDL::getId("Lambda_b0");
 382
 383   // Anti-Baryons
 384   static EvtId LAMCM=EvtPDL::getId("anti-Lambda_c-");
 385   static EvtId LAMC1M=EvtPDL::getId("anti-Lambda_c(2593)-");
 386   static EvtId LAMC2M=EvtPDL::getId("anti-Lambda_c(2625)-");
 387   static EvtId LAMBB=EvtPDL::getId("anti-Lambda_b0");
 388
 389   // Set the spin density matrix of the parent baryon
 390   EvtSpinDensity rho;
 391   rho.setDim(2);
 392   rho.set(0,0,r00);
 393   rho.set(0,1,r01);
 394
 395   rho.set(1,0,r10);
 396   rho.set(1,1,r11);
 397
 398   EvtVector4R vector4P = parent->getP4Lab();
 399   double pmag = vector4P.d3mag();
 400   double cosTheta = vector4P.get(3)/pmag;
 401
 402   double theta = acos(cosTheta);
 403   double phi = atan2(vector4P.get(2), vector4P.get(1));
 404
 405   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   EvtVector4R p4b;
 410   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   EvtVector4R p4daught = parent->getDaug(0)->getP4();
 414
 415   // Add the lepton and neutrino 4 momenta to find q (q^2)
 416   EvtVector4R q = parent->getDaug(1)->getP4()
 417     + parent->getDaug(2)->getP4();
 418
 419   double q2 = q.mass2();
 420
 421
 422   EvtId l_num = parent->getDaug(1)->getId();
 423   EvtId bar_num = parent->getDaug(0)->getId();
 424   EvtId par_num = parent->getId();
 425
 426   double baryonmass = parent->getDaug(0)->mass();
 427
 428   // Handle spin-1/2 daughter baryon Dirac spinor cases
 429   if( EvtPDL::getSpinType(parent->getDaug(0)->getId())==EvtSpinType::DIRAC ) {
 430
 431     // Set the form factors
 432     double f1,f2,f3,g1,g2,g3;
 433     FormFactors->getdiracff( par_num,
 434                              bar_num,
 435                              q2,
 436                              baryonmass,
 437                              &f1, &f2, &f3,
 438                              &g1, &g2, &g3);
 439
 440     const double form_fact[6] = {f1, f2, f3, g1, g2, g3};
 441
 442     EvtVector4C b11, b12, b21, b22, l1, l2;
 443
 444     //  Lepton Current
 445     if(l_num==EM || l_num==MUM || l_num==TAUM){
 446
 447       l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
 448                             parent->getDaug(2)->spParentNeutrino());
 449       l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
 450                             parent->getDaug(2)->spParentNeutrino());
 451
 452     } else if (l_num==EP || l_num==MUP || l_num==TAUP) {
 453
 454       l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
 455                             parent->getDaug(1)->spParent(0));
 456       l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
 457                             parent->getDaug(1)->spParent(1));
 458
 459     } else {
 460       report(ERROR,"EvtGen")<< "Wrong lepton number \n";
 461       ::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     if ( (par_num==LAMB && bar_num==LAMCP)
 476          || (par_num==LAMBB && bar_num==LAMCM) ) {
 477
 478       // Set particle/anti-particle flag
 479       if (bar_num==LAMCP)
 480         pflag = 0;
 481       else if (bar_num==LAMCM)
 482         pflag = 2;
 483
 484       b11=EvtBaryonVACurrent(parent->getDaug(0)->spParent(0),
 485                              parent->sp(0),
 486                              p4b, p4daught, form_fact, pflag);
 487       b21=EvtBaryonVACurrent(parent->getDaug(0)->spParent(0),
 488                              parent->sp(1),
 489                              p4b, p4daught, form_fact, pflag);
 490       b12=EvtBaryonVACurrent(parent->getDaug(0)->spParent(1),
 491                              parent->sp(0),
 492                              p4b, p4daught, form_fact, pflag);
 493       b22=EvtBaryonVACurrent(parent->getDaug(0)->spParent(1),
 494                              parent->sp(1),
 495                              p4b, p4daught, form_fact, pflag);
 496     }
 497
 498     // Handle 1/2+ -> 1/2- second
 499     else if( (par_num==LAMB && bar_num==LAMC1P)
 500              || (par_num==LAMBB && bar_num==LAMC1M) ) {
 501
 502       // Set particle/anti-particle flag
 503       if (bar_num==LAMC1P)
 504         pflag = 1;
 505       else if (bar_num==LAMC1M)
 506         pflag = 3;
 507
 508       b11=EvtBaryonVACurrent((parent->getDaug(0)->spParent(0)),
 509                              (EvtGammaMatrix::g5()*parent->sp(0)),
 510                              p4b, p4daught, form_fact, pflag);
 511       b21=EvtBaryonVACurrent((parent->getDaug(0)->spParent(0)),
 512                              (EvtGammaMatrix::g5()*parent->sp(1)),
 513                              p4b, p4daught, form_fact, pflag);
 514       b12=EvtBaryonVACurrent((parent->getDaug(0)->spParent(1)),
 515                              (EvtGammaMatrix::g5()*parent->sp(0)),
 516                              p4b, p4daught, form_fact, pflag);
 517       b22=EvtBaryonVACurrent((parent->getDaug(0)->spParent(1)),
 518                              (EvtGammaMatrix::g5()*parent->sp(1)),
 519                              p4b, p4daught, form_fact, pflag);
 520
 521     }
 522
 523     else {
 524       report(ERROR,"EvtGen") << "Dirac semilep. baryon current "
 525                              << "not implemented for this decay sequence."
 526                              << std::endl;
 527       ::abort();
 528     }
 529
 530     amp.vertex(0,0,0,l1*b11);
 531     amp.vertex(0,0,1,l2*b11);
 532
 533     amp.vertex(1,0,0,l1*b21);
 534     amp.vertex(1,0,1,l2*b21);
 535
 536     amp.vertex(0,1,0,l1*b12);
 537     amp.vertex(0,1,1,l2*b12);
 538
 539     amp.vertex(1,1,0,l1*b22);
 540     amp.vertex(1,1,1,l2*b22);
 541
 542   }
 543
 544   // Need special handling for the spin-3/2 daughter baryon
 545   // Rarita-Schwinger spinor cases
 546   else if( EvtPDL::getSpinType(parent->getDaug(0)->getId())==EvtSpinType::RARITASCHWINGER ) {
 547
 548     // Set the form factors
 549     double f1,f2,f3,f4,g1,g2,g3,g4;
 550     FormFactors->getraritaff( par_num,
 551                               bar_num,
 552                               q2,
 553                               baryonmass,
 554                               &f1, &f2, &f3, &f4,
 555                               &g1, &g2, &g3, &g4);
 556
 557     const double form_fact[8] = {f1, f2, f3, f4, g1, g2, g3, g4};
 558
 559     EvtId l_num = parent->getDaug(1)->getId();
 560
 561     EvtVector4C b11, b12, b21, b22, b13, b23, b14, b24, l1, l2;
 562
 563     //  Lepton Current
 564     if (l_num==EM || l_num==MUM || l_num==TAUM) {
 565       //  Lepton Current
 566       l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
 567                             parent->getDaug(2)->spParentNeutrino());
 568       l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
 569                             parent->getDaug(2)->spParentNeutrino());
 570     }
 571     else if (l_num==EP || l_num==MUP || l_num==TAUP) {
 572       l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
 573                             parent->getDaug(1)->spParent(0));
 574       l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
 575                             parent->getDaug(1)->spParent(1));
 576
 577     } else {
 578       report(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     if (par_num==LAMB && bar_num==LAMC2P) {
 589       // Set flag for particle case
 590       pflag = 0;
 591     }
 592     else if (par_num==LAMBB && bar_num==LAMC2M) {
 593       // Set flag for anti-particle case
 594       pflag = 1;
 595     }
 596     else {
 597       report(ERROR,"EvtGen") << "Rarita-Schwinger semilep. baryon current "
 598                              << "not implemented for this decay sequence."
 599                              << std::endl;
 600       ::abort();
 601     }
 602
 603     // Baryon current
 604     b11=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(0),
 605                                  parent->sp(0),
 606                                  p4b, p4daught, form_fact, pflag);
 607     b21=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(0),
 608                                  parent->sp(1),
 609                                  p4b, p4daught, form_fact, pflag);
 610
 611     b12=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(1),
 612                                  parent->sp(0),
 613                                  p4b, p4daught, form_fact, pflag);
 614     b22=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(1),
 615                                  parent->sp(1),
 616                                  p4b, p4daught, form_fact, pflag);
 617
 618     b13=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(2),
 619                                  parent->sp(0),
 620                                  p4b, p4daught, form_fact, pflag);
 621     b23=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(2),
 622                                  parent->sp(1),
 623                                  p4b, p4daught, form_fact, pflag);
 624
 625     b14=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(3),
 626                                  parent->sp(0),
 627                                  p4b, p4daught, form_fact, pflag);
 628     b24=EvtBaryonVARaritaCurrent(parent->getDaug(0)->spRSParent(3),
 629                                  parent->sp(1),
 630                                  p4b, p4daught, form_fact, pflag);
 631
 632     amp.vertex(0,0,0,l1*b11);
 633     amp.vertex(0,0,1,l2*b11);
 634
 635     amp.vertex(1,0,0,l1*b21);
 636     amp.vertex(1,0,1,l2*b21);
 637
 638     amp.vertex(0,1,0,l1*b12);
 639     amp.vertex(0,1,1,l2*b12);
 640
 641     amp.vertex(1,1,0,l1*b22);
 642     amp.vertex(1,1,1,l2*b22);
 643
 644     amp.vertex(0,2,0,l1*b13);
 645     amp.vertex(0,2,1,l2*b13);
 646
 647     amp.vertex(1,2,0,l1*b23);
 648     amp.vertex(1,2,1,l2*b23);
 649
 650     amp.vertex(0,3,0,l1*b14);
 651     amp.vertex(0,3,1,l2*b14);
 652
 653     amp.vertex(1,3,0,l1*b24);
 654     amp.vertex(1,3,1,l2*b24);
 655
 656   }
 657
 658 }
 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   EvtComplex cv = EvtComplex(1.0, 0.);
 675   EvtComplex ca = EvtComplex(1.0, 0.);
 676
 677   EvtComplex cg0 = EvtComplex(1.0, 0.);
 678   EvtComplex cg5 = EvtComplex(1.0, 0.);
 679
 680   // antiparticle- same parity parent & daughter
 681   if( pflag == 2 ) {
 682     cv = EvtComplex(-1.0, 0.);
 683     ca = EvtComplex(1.0, 0.);
 684
 685     cg0 =  EvtComplex(1.0, 0.0);
 686     cg5 =  EvtComplex(0.0, -1.0);
 687   }
 688   // antiparticle- opposite parity parent & daughter
 689   else if( pflag == 3) {
 690     cv = EvtComplex(1.0, 0.);
 691     ca = EvtComplex(-1.0, 0.);
 692
 693     cg0 =  EvtComplex(0.0, -1.0);
 694     cg5 =  EvtComplex(1.0, 0.0);
 695   }
 696
 697   EvtVector4C t[6];
 698
 699   // Term 1 = \bar{u}(p',s')*(F_1(q^2)*\gamma_{mu})*u(p,s)
 700   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   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   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   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   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   t[5] = cg5*EvtLeptonPCurrent( Bf, Bi ) * (daught/daught.mass());
 716
 717   // Sum the individual terms
 718   EvtVector4C current = (ff[0]*t[0] + ff[1]*t[1] + ff[2]*t[2]
 719                          - ff[3]*t[3] - ff[4]*t[4] - ff[5]*t[5]);
 720
 721   return current;
 722 }
 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   EvtComplex cv = EvtComplex(1.0, 0.);
 736   EvtComplex ca = EvtComplex(1.0, 0.);
 737
 738   EvtComplex cg0 = EvtComplex(1.0, 0.);
 739   EvtComplex cg5 = EvtComplex(1.0, 0.);
 740
 741   // antiparticle
 742   if( pflag == 1 ) {
 743     cv = EvtComplex(-1.0, 0.);
 744     ca = EvtComplex(1.0, 0.);
 745
 746     cg0 =  EvtComplex(1.0, 0.0);
 747     cg5 =  EvtComplex(0.0, -1.0);
 748  }
 749
 750   EvtVector4C t[8];
 751   EvtTensor4C id;
 752   id.setdiag(1.0,1.0,1.0,1.0);
 753
 754   EvtDiracSpinor tmp;
 755   for(int i=0;i<4;i++){
 756     tmp.set_spinor(i,Bf.getVector(i)*parent);
 757   }
 758
 759   EvtVector4C v1,v2;
 760   for(int i=0;i<4;i++){
 761     v1.set(i,EvtLeptonSCurrent(Bf.getSpinor(i),Bi));
 762     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   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   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   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   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   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   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   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   t[7] = cg5*(id.cont2(v2));
 795
 796   // Sum the individual terms
 797   EvtVector4C current = (ff[0]*t[0] + ff[1]*t[1] + ff[2]*t[2] + ff[3]*t[3]
 798                          - ff[4]*t[4] - ff[5]*t[5] - ff[6]*t[6] - ff[7]*t[7]);
 799
 800   return current;
 801 }