TEvtGen/EvtGenModels/EvtLb2Lll.cpp

   1 //----------------------------------------------------------------------------------
   2 //
   3 // Module: EvtLb2Lll.cpp
   4 //
   5 // Desription: Routine to implement Lambda_b0 -> Lambda_0 l+ l- decays accroding to
   6 //             several models: Chen. Geng.
   7 //                             Aliev. Ozpineci. Savci.
   8 //
   9 // Modification history:
  10 //
  11 //  10/07/2012  MK   Fix calculation of N1, N2; based on hep-ph/021144
  12 //  09/02/2009  PR   Commented check for (anti-)Lambda0 names
  13 //  15/09/2004  PR   Module created according to PHSP model
  14 //  20/02/2005  PR   Added parameters, created matrix element (without polarization)
  15 //  04/03/2005  PR   LD contrib., corrected WC eff. according to Chen. Geng.
  16 //
  17 // Todo list:
  18 //
  19 //   - Properly handle antiparticles, needs change of u, ubar to v, vbar in
  20 //     hadronic current, or other way of putting that in
  21 //----------------------------------------------------------------------------------
  22
  23 #ifdef WIN32
  24 #pragma warning( disable : 4786 )
  25 // Disable anoying warning about symbol size
  26 #endif
  27
  28 #include "EvtGenModels/EvtLb2Lll.hh"
  29 #include "EvtGenBase/EvtParticle.hh"
  30 #include "EvtGenBase/EvtDiracParticle.hh"
  31 #include "EvtGenBase/EvtPDL.hh"
  32 #include "EvtGenBase/EvtDiracSpinor.hh"
  33 #include "EvtGenBase/EvtTensor4C.hh"
  34 #include "EvtGenBase/EvtVector4C.hh"
  35 #include "EvtGenBase/EvtVector4R.hh"
  36 #include "EvtGenBase/EvtComplex.hh"
  37 #include "EvtGenBase/EvtGammaMatrix.hh"
  38 #include <stdio.h>
  39 #include <string.h>
  40
  41 EvtLb2Lll::~EvtLb2Lll() {}
  42
  43 EvtDecayBase* EvtLb2Lll::clone(){
  44   return new EvtLb2Lll;
  45 }
  46
  47 std::string EvtLb2Lll::getName(){
  48   return "Lb2Lll";
  49 }
  50
  51 void EvtLb2Lll::init(){
  52
  53   if(getNArg()>8){ // Decay parameters
  54     report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected max. 8 arguments but found: " << getNArg() << std::endl;
  55     report(INFO ,"EvtGen") << "  1. Lambda_b0 polarization - zero is default" << std::endl;
  56     report(INFO ,"EvtGen") << "  2. Model type - \"SM\" for Standard Model is default" << std::endl;
  57     report(INFO ,"EvtGen") << "  3. Form-Factors - \"HQET\" is used by default" << std::endl;
  58     report(INFO ,"EvtGen") << "  4. How to set polarization - \"ModifiedSpinors\" is default" << std::endl;
  59     report(INFO ,"EvtGen") << "  5. Include long distance (LD) effects - \"SD\" (no) is default" << std::endl;
  60     report(INFO ,"EvtGen") << "  6. NonFactorizable contribution (omega) to b->sg decay at q2=0 " << std::endl;
  61     report(INFO ,"EvtGen") << "  7. Note on every x-th decay" << std::endl;
  62     report(INFO ,"EvtGen") << "  8. Maximum probability - automatic by default" << std::endl;
  63     ::abort();
  64   }
  65
  66   if(getNDaug()!=3){ // Check that there are 3 daughters only
  67     report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected 3 daughters but found: " << getNDaug() << std::endl;
  68     ::abort();
  69   }
  70
  71   // TODO: better check based on spin and falvour is needed to allow usage of aliases !
  72   if(EvtPDL::name(getParentId())=="Lambda_b0"){ // Check daughters of Lambda_b0
  73     report(INFO,"EvtGen") << " EvtLb2Lll generator found Lambda_b0" << std::endl;
  74     //if(EvtPDL::name(getDaug(0))!="Lambda0"){
  75     //  report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected Lambda0 daughter but found: " << EvtPDL::name(getDaug(0)) << std::endl;
  76     //  ::abort();
  77     //}
  78     if(EvtPDL::name(getDaug(1))=="e-" && EvtPDL::name(getDaug(2))=="e+"){
  79       m_decayName="Lambda_b0 -> Lambda0 e- e+";
  80       report(INFO,"EvtGen") << " EvtLb2Lll generator found decay:  Lambda_b0 -> Lambda0 e- e+" << std::endl;
  81     }else if(EvtPDL::name(getDaug(1))=="mu-" && EvtPDL::name(getDaug(2))=="mu+"){
  82       m_decayName="Lambda_b0 -> Lambda0 mu- mu+";
  83       report(INFO,"EvtGen") << " EvtLb2Lll generator found decay:  Lambda_b0 -> Lambda0 mu- mu+" << std::endl;
  84     }else if(EvtPDL::name(getDaug(1))=="tau-" && EvtPDL::name(getDaug(2))=="tau+"){
  85       m_decayName="Lambda_b0 -> Lambda0 tau- tau+";
  86       report(INFO,"EvtGen") << " EvtLb2Lll generator found decay:  Lambda_b0 -> Lambda0 tau- tau+" << std::endl;
  87     }else{
  88       report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected lepton pair daughters but found:  " << EvtPDL::name(getDaug(1)) << " " << EvtPDL::name(getDaug(2)) << std::endl;
  89       ::abort();
  90     }
  91   //TODO: The model is known not to work correctly for anti-Lambda_b0 (A_FB does not change its sign)
  92   }else if(EvtPDL::name(getParentId())=="anti-Lambda_b0"){ // Check daughters of anti-Lambda_b0
  93     report(INFO,"EvtGen") << " EvtLb2Lll generator found anti-Lambda_b0" << std::endl;
  94     //if(EvtPDL::name(getDaug(0))!="anti-Lambda0"){
  95     //  report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected anti-Lambda0 daughter but found: " << EvtPDL::name(getDaug(0)) << std::endl;
  96     //  ::abort();
  97     //}
  98     if(EvtPDL::name(getDaug(1))=="e+" && EvtPDL::name(getDaug(2))=="e-"){
  99       m_decayName="anti-Lambda_b0 -> anti-Lambda0 e+ e-";
 100       report(INFO,"EvtGen") << " EvtLb2Lll generator found decay:  anti-Lambda_b0 -> anti-Lambda0 e+ e-" << std::endl;
 101     }else if(EvtPDL::name(getDaug(1))=="mu+" && EvtPDL::name(getDaug(2))=="mu-"){
 102       m_decayName="anti-Lambda_b0 -> anti-Lambda0 mu+ mu-";
 103       report(INFO,"EvtGen") << " EvtLb2Lll generator found decay:  anti-Lambda_b0 -> anti-Lambda0 mu+ mu-" << std::endl;
 104     }else if(EvtPDL::name(getDaug(1))=="tau-" && EvtPDL::name(getDaug(2))=="tau+"){
 105       m_decayName="anti-Lambda_b0 -> anti-Lambda0 tau+ tau-";
 106       report(INFO,"EvtGen") << " EvtLb2Lll generator found decay:  anti-Lambda_b0 -> anti-Lambda0 tau+ tau-" << std::endl;
 107     }else{
 108       report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected lepton pair daughters but found:  " << EvtPDL::name(getDaug(1)) << " " << EvtPDL::name(getDaug(2)) << std::endl;
 109       ::abort();
 110     }
 111   }else{ // This model is not intended for decay of anything else than (anti-)Lambda_b0
 112     report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected (anti-)Lambda_b0 parent but found:  " << EvtPDL::name(getParentId()) << std::endl;
 113     ::abort();
 114   }
 115
 116   // Read and check all parameters
 117   if(getNArg()>0){
 118     if(getArg(0)>1. || getArg(0)<-1.){
 119       report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll expects polarization to be in interval <-1,1>, not " << getArg(0) << std::endl;
 120       ::abort();
 121     }
 122     m_polarizationLambdab0 = getArg(0);
 123   }else{
 124     m_polarizationLambdab0 = 0;
 125   }
 126   report(INFO,"EvtGen") << " EvtLb2Lll set Lambda_b0 polarization to " << m_polarizationLambdab0 << std::endl;
 127
 128   if(getNArg()>1){
 129     if(getArgStr(1).substr(1,getArgStr(1).size()-2)!="SM" &&
 130        getArgStr(1).substr(1,getArgStr(1).size()-2)!="-C7_SM" &&
 131        getArgStr(1).substr(1,getArgStr(1).size()-2)!="SUSY-ChenGeng"){
 132       report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll doesn't know this physics model: " << getArgStr(1) << std::endl;
 133       ::abort();
 134     }
 135     m_HEPmodel = getArgStr(1).substr(1,getArgStr(1).size()-2);
 136   }else{
 137     m_HEPmodel = "SM";
 138   }
 139   report(INFO,"EvtGen") << " EvtLb2Lll will use this physics model: " << m_HEPmodel << std::endl;
 140
 141   if(getNArg()>2){
 142     if(getArgStr(2).substr(1,getArgStr(2).size()-2)!="HQET" &&
 143        getArgStr(2).substr(1,getArgStr(2).size()-2)!="HQET-noF2" &&
 144        getArgStr(2).substr(1,11)                   !="HQET-delta="){
 145       report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll doesn't know this Form-Factors model: " << getArgStr(2) << std::endl;
 146       ::abort();
 147     }
 148     m_FFtype = getArgStr(2).substr(1,getArgStr(2).size()-2);
 149   }else{
 150     m_FFtype = "HQET";
 151   }
 152   report(INFO,"EvtGen") << " EvtLb2Lll will use this Form-Factors model: " << m_FFtype << std::endl;
 153
 154   if(getNArg()>3){
 155     if(getArgStr(3).substr(1,getArgStr(3).size()-2)!="Unpolarized"){
 156       report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll doesn't know kind of introducing polarization: " << getArgStr(3) << std::endl;
 157       ::abort();
 158     }
 159     m_polarizationIntroduction = getArgStr(3).substr(1,getArgStr(3).size()-2);
 160   }else{
 161     m_polarizationIntroduction = "Unpolarized";
 162   }
 163   report(INFO,"EvtGen") << " EvtLb2Lll will use this kind of introducing polarization: " << m_polarizationIntroduction << std::endl;
 164
 165   if(getNArg()>4){
 166     if(getArgStr(4).substr(1,getArgStr(4).size()-2)!="SD" && getArgStr(4).substr(1,getArgStr(4).size()-2)!="LD"){
 167       report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll didn't find SD or LD parameter: " << getArgStr(4) << std::endl;
 168       ::abort();
 169     }
 170     m_effectContribution = getArgStr(5).substr(1,getArgStr(4).size()-2);
 171   }else{
 172     m_effectContribution = "SD";
 173   }
 174   report(INFO,"EvtGen") << " EvtLb2Lll will include contribution from these effects: " << m_effectContribution << std::endl;
 175
 176   if(getNArg()>5){
 177     if(fabs(getArg(5))>0.15){
 178       report(WARNING,"EvtGen") << " WARNING: EvtLb2Lll found very high contribution to b->sg decay at q2=0: " << getArg(5) << std::endl;
 179     }
 180     m_omega = getArg(5);
 181   }else{
 182     m_omega = 0;
 183   }
 184   report(INFO,"EvtGen") << " EvtLb2Lll will use this contribution to b->sg decay at q2=0: " << m_omega << std::endl;
 185
 186   if(getNArg()>6) m_noTries = (long)(getArg(6));
 187   else            m_noTries = 0;
 188
 189   if(getNArg()>7){
 190     if(getArg(7)<0.){
 191       report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll expects positive maximum probability not : " << getArg(7) << std::endl;
 192       ::abort();
 193     }
 194     m_maxProbability = getArg(7);
 195   }else{
 196     m_maxProbability = 0.;
 197   }
 198   report(INFO,"EvtGen") << " EvtLb2Lll maximum probability was set to " << m_maxProbability << std::endl;
 199   m_poleSize=0;
 200
 201   // Initialize Wilson coefficients by Buras and Munz
 202   // TODO: should have common W.C. source for all decays in EvtGen
 203   m_WC.CalculateAllCoefficients();
 204
 205 }
 206
 207 void EvtLb2Lll::initProbMax(){
 208
 209   report(INFO,"EvtGen") << " EvtLb2Lll is finding maximum probability ... " << std::endl;
 210
 211   if(m_maxProbability==0){
 212
 213     EvtDiracParticle *parent = new EvtDiracParticle;
 214     parent->noLifeTime();
 215     parent->init(getParentId(),EvtVector4R(EvtPDL::getMass(getParentId()),0,0,0));
 216     parent->setDiagonalSpinDensity();
 217
 218     EvtAmp amp;
 219     EvtId daughters[3] = {getDaug(0),getDaug(1),getDaug(2)};
 220     amp.init(getParentId(),3,daughters);
 221     parent->makeDaughters(3,daughters);
 222     EvtParticle *lambda = parent->getDaug(0);
 223     EvtParticle *lep1   = parent->getDaug(1);
 224     EvtParticle *lep2   = parent->getDaug(2);
 225     lambda -> noLifeTime();
 226     lep1   -> noLifeTime();
 227     lep2   -> noLifeTime();
 228
 229     EvtSpinDensity rho;
 230     rho.setDiag(parent->getSpinStates());
 231
 232     double M0 = EvtPDL::getMass(getParentId());
 233     double mL = EvtPDL::getMass(getDaug(0));
 234     double m1 = EvtPDL::getMass(getDaug(1));
 235     double m2 = EvtPDL::getMass(getDaug(2));
 236
 237     double q2,pstar,elambda,theta;
 238     double q2min = (m1+m2)*(m1+m2);
 239     double q2max = (M0-mL)*(M0-mL);
 240
 241     EvtVector4R p4lambda,p4lep1,p4lep2,boost;
 242
 243     report(INFO,"EvtGen") << " EvtLb2Lll is probing whole phase space ..." << std::endl;
 244
 245     int i,j;
 246     double prob=0;
 247     for(i=0;i<=100;i++){
 248       q2 = q2min+i*(q2max-q2min)/100.;
 249       elambda = (M0*M0+mL*mL-q2)/2/M0;
 250       if(i==0) pstar = 0;
 251       else     pstar = sqrt(q2-(m1+m2)*(m1+m2))*sqrt(q2-(m1-m2)*(m1-m2))/2/sqrt(q2);
 252       boost.set(M0-elambda,0,0,+sqrt(elambda*elambda-mL*mL));
 253       p4lambda.set(elambda,0,0,-sqrt(elambda*elambda-mL*mL));
 254       for(j=0;j<=45;j++){
 255         theta = j*EvtConst::pi/45;
 256         p4lep1.set(sqrt(pstar*pstar+m1*m1),0,+pstar*sin(theta),+pstar*cos(theta));
 257         p4lep2.set(sqrt(pstar*pstar+m2*m2),0,-pstar*sin(theta),-pstar*cos(theta));
 258         //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " pstar: " << pstar << std::endl;
 259         p4lep1 = boostTo(p4lep1,boost);
 260         p4lep2 = boostTo(p4lep1,boost);
 261         lambda -> init(getDaug(0),p4lambda);
 262         lep1   -> init(getDaug(1),p4lep1  );
 263         lep2   -> init(getDaug(2),p4lep2  );
 264         calcAmp(&amp,parent);
 265         prob = rho.normalizedProb(amp.getSpinDensity());
 266         //std::cout << "q2:  " << q2 << " \t theta:  " << theta << " \t prob:  " << prob << std::endl;
 267         //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " q2-q2min: " << q2-(m1+m2)*(m1+m2) << std::endl;
 268         if(prob>m_maxProbability){
 269           report(INFO,"EvtGen") << "  - probability " << prob << " found at q2 = " << q2 << " (" << 100*(q2-q2min)/(q2max-q2min) << " %) and theta = " << theta*180/EvtConst::pi << std::endl;
 270           m_maxProbability=prob;
 271         }
 272       }
 273       //::abort();
 274     }
 275
 276     //m_poleSize = 0.04*q2min;
 277     m_maxProbability *= 1.2;
 278     delete parent;
 279   }
 280
 281   setProbMax(m_maxProbability);
 282   report(INFO,"EvtGen") << " EvtLb2Lll set up maximum probability to " << m_maxProbability << std::endl;
 283
 284 }
 285
 286 void EvtLb2Lll::decay(EvtParticle* parent){
 287
 288   //setWeight(parent->initializePhaseSpace(getNDaug(),getDaugs(),m_poleSize,1,2));
 289   parent->initializePhaseSpace(getNDaug(),getDaugs());
 290   calcAmp(&_amp2,parent);
 291
 292 }
 293
 294 void EvtLb2Lll::calcAmp(EvtAmp *amp,EvtParticle *parent){
 295
 296   static long noTries=0;
 297   static double delta=0;
 298
 299   EvtComplex Matrix[2][2][2][2];
 300
 301   EvtComplex i1(0,1);
 302
 303   int i,j,spins[4];
 304   char ch;
 305
 306   double r,M_L,M_Lb,M_s,M_c,M_b,q2,alpha,M_W,M_t;
 307   double M_psi[2]={0,0},Gamma_psi[2]={0,0},k_psi[2]={0,0};
 308   double F0_1,F0_2,a_F1,a_F2,b_F1,b_F2,F1,F2;
 309   double f_1,f_2,f_3,g_1,g_2,g_3,f_1T,f_2T,f_3T,g_1T,g_2T,g_3T,f_TV,f_TS,g_TV(0.0),g_TS,f_T,g_T;
 310   EvtComplex A1,A2,A3,B1,B2,B3,D1,D2,D3,E1,E2,E3,N1,N2,H1,H2;
 311   EvtComplex C_SL,C_BR,C_LLtot,C_LRtot,C_LL,C_LR,C_RL,C_RR,C_LRLR,C_RLLR,C_LRRL,C_RLRL,C_T,C_TE;
 312   EvtComplex Yld,C_7eff,C_9eff;
 313   EvtComplex V_ts,V_tb;
 314
 315   EvtVector4C lbar_Gmu_l[2][2],lbar_GmuG5_l[2][2],hbar_GmuPlusG5_h[2][2],hbar_GmuMinusG5_h[2][2],hbar_Gmu_h[2][2];
 316   EvtComplex  lbar_l[2][2],lbar_G5_l[2][2],hbar_1PlusG5_h[2][2],hbar_1MinusG5_h[2][2],hbar_G5_h[2][2],hbar_h[2][2];
 317   EvtTensor4C lbar_Smunu_l[2][2],lbar_ESmunu_l[2][2],hbar_SmunuPlusG5_h[2][2],hbar_SmunuMinusG5_h[2][2],hbar_Smunu_h[2][2];
 318   EvtVector4R q_mu,P_mu;
 319
 320   EvtDiracSpinor parent__spParent[2];
 321
 322   M_Lb    = parent->mass();
 323   M_L     = parent->getDaug(0)->mass();
 324   M_s     = 0.13;
 325   M_c     = 1.35;
 326   M_b     = 4.8;
 327   alpha   = 1./137.036;
 328   M_W     = 80.425;
 329   M_t     = 174.3;
 330   M_psi[0] = 3.096916;
 331   M_psi[1] = 3.686093;
 332   if(m_decayName=="Lambda_b0 -> Lambda0 e- e+" || m_decayName=="anti-Lambda_b0 -> anti-Lambda0 e+ e-"){
 333     Gamma_psi[0] = 5.40;
 334     Gamma_psi[1] = 2.12;
 335   }
 336   if(m_decayName=="Lambda_b0 -> Lambda0 mu- mu+" || m_decayName=="anti-Lambda_b0 -> anti-Lambda0 mu+ mu-"){
 337     Gamma_psi[0] = 5.35;
 338     Gamma_psi[1] = 2.05;
 339   }
 340   if(m_decayName=="Lambda_b0 -> Lambda0 tau- tau+" || m_decayName=="anti-Lambda_b0 -> anti-Lambda0 tau+ tau-"){
 341     Gamma_psi[0] = 0.00;
 342     Gamma_psi[1] = 0.79;
 343   }
 344   if(m_effectContribution=="LD"){
 345     k_psi[0] = 1.65;
 346     k_psi[1] = 1.65;
 347   }
 348   //G_F   = 1.16637e-5;
 349   //V_tb  = sqrt(1-pow(0.0413,2))*sqrt(1-pow(0.0037,2));
 350   //V_ts  = -sqrt(1-pow(0.2243,2))*0.0413-0.2243*sqrt(1-pow(0.0413,2))*0.0037*(cos(60*EvtConst::pi/180)+i1*sin(60*EvtConst::pi/180));
 351
 352   P_mu = parent->getP4Restframe()+parent->getDaug(0)->getP4();
 353   q_mu = parent->getP4Restframe()-parent->getDaug(0)->getP4();
 354   q2   = q_mu.mass2();
 355
 356   if(m_noTries>0) if(!((++noTries)%m_noTries)) report(DEBUG,"EvtGen") << " EvtLb2Lll already finished " << noTries << " matrix element calculations" << std::endl;
 357
 358   if(m_FFtype=="HQET"){
 359     r = M_L*M_L/M_Lb/M_Lb;
 360     F0_1 = +0.462;
 361     F0_2 = -0.077;
 362     a_F1 = -0.0182;
 363     a_F2 = -0.0685;
 364     b_F1 = -0.000176;
 365     b_F2 = +0.001460;
 366     F1 = F0_1/(1.0-(q2/M_Lb/M_Lb)*(a_F1-b_F1*(q2/M_Lb/M_Lb)));
 367     F2 = F0_2/(1.0-(q2/M_Lb/M_Lb)*(a_F2-b_F2*(q2/M_Lb/M_Lb)));
 368     g_1 = f_1 = f_2T = g_2T = F1+sqrt(r)*F2;
 369     //std::cout << " F1: " << F1 << "  F2: " << F2 << "  r: " << r << "  M_L: " << M_L << "  M_Lb: " << M_Lb << std::endl;
 370     //std::cout << " sqrt(q2): " << sqrt(q2) << "  q2: " << q2 << "  M_Lb^2" << M_Lb*M_Lb << std::endl;
 371     g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = F2/M_Lb;
 372     g_TS = f_TS = 0;
 373     g_1T = f_1T = F2/M_Lb*q2;
 374     g_3T = +F2/M_Lb*(M_Lb+M_L);
 375     f_3T = -F2/M_Lb*(M_Lb-M_L);
 376     f_T = f_2T-f_TS*q2;
 377     g_T = g_2T-g_TS*q2;
 378   }else if(strstr(m_FFtype.c_str(),"HQET-delta=")==m_FFtype.c_str()){
 379     //report(WARNING,"EvtGen") << " WARNING: HQET-delta FF model should be checked for correctness" << std::endl;
 380     if(delta==0) sscanf(m_FFtype.c_str(),"%c%c%c%c%c%c%c%c%c%c%c%lf",&ch,&ch,&ch,&ch,&ch,&ch,&ch,&ch,&ch,&ch,&ch,&delta);
 381     r = M_L*M_L/M_Lb/M_Lb;
 382     F0_1 = +0.462;
 383     F0_2 = -0.077;
 384     a_F1 = -0.0182;
 385     a_F2 = -0.0685;
 386     b_F1 = -0.000176;
 387     b_F2 = +0.001460;
 388     F1 = F0_1/(1.0-(q2/M_Lb/M_Lb)*(a_F1-b_F1*(q2/M_Lb/M_Lb)));
 389     F2 = F0_2/(1.0-(q2/M_Lb/M_Lb)*(a_F2-b_F2*(q2/M_Lb/M_Lb)));
 390     g_1 = f_1 = f_2T = g_2T = F1+sqrt(r)*F2;
 391     g_1 += delta*g_1;
 392     f_1 -= delta*f_1;
 393     g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = F2/M_Lb;
 394     g_TS = f_TS = 0;
 395     g_1T = f_1T = F2/M_Lb*q2;
 396     g_3T = +F2/M_Lb*(M_Lb+M_L);
 397     f_3T = -F2/M_Lb*(M_Lb-M_L);
 398     f_T = f_2T-f_TS*q2;
 399     g_T = g_2T-g_TS*q2;
 400   }else if(m_FFtype=="HQET-noF2"){
 401     //report(WARNING,"EvtGen") << " WARNING: HQET-noF2 FF model should be checked for correctness" << std::endl;
 402     r = M_L*M_L/M_Lb/M_Lb;
 403     F0_1 = +0.462;
 404     a_F1 = -0.0182;
 405     b_F1 = -0.000176;
 406     F1 = F0_1/(1.0-(q2/M_Lb/M_Lb)*(a_F1-b_F1*(q2/M_Lb/M_Lb)));
 407     g_1 = f_1 = f_2T = g_2T = F1;
 408     g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = 0;
 409     g_TS = f_TS = 0;
 410     g_1T = f_1T = 0;
 411     g_3T = 0;
 412     f_3T = 0;
 413     f_T = f_2T-f_TS*q2;
 414     g_T = g_2T-g_TS*q2;
 415   }else{ // general relations for Form-Factors
 416     f_1 = f_2 = f_3 = g_1 = g_2 = g_3 = f_3T = g_3T = f_TS = g_TS = f_T = g_T = f_TV = 0;
 417     f_2T = f_T+f_TS*q2;
 418     f_1T = (f_TV+f_TS*(M_L+M_Lb))*q2;
 419     f_1T = -q2/(M_Lb-M_L)*f_3T;
 420     g_2T = g_T+g_TS*q2;
 421     g_1T = (g_TV-g_TS*(M_L-M_Lb))*q2;
 422     g_1T = +q2/(M_Lb+M_L)*g_3T;
 423     report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll - unknown Form-Factors model: " << m_FFtype << " - this should never happen !" << std::endl;
 424     ::abort();
 425   }
 426
 427   if(m_HEPmodel=="SM"){
 428     C_LL=C_LR=C_RL=C_RR=C_LRLR=C_RLLR=C_LRRL=C_RLRL=C_T=C_TE=EvtComplex(0,0);
 429     Yld = m_WC.Yld(q2,k_psi,Gamma_psi,M_psi,2,m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),1./alpha);
 430     C_7eff = m_WC.GetC7eff0() + m_WC.C7b2sg(m_WC.GetStrongCouplingConst(),m_WC.GetEta(),m_WC.GetC2(),M_t,M_W) + m_omega*(m_WC.hzs(M_c/M_b,q2/M_b/M_b,M_b,M_b) + Yld);
 431     C_9eff = Yld + m_WC.C9efftilda(M_c/M_b,q2/M_b/M_b,m_WC.GetStrongCouplingConst(),m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),m_WC.GetC9tilda(),m_WC.GetRenormSchemePar());
 432     C_SL = -2*M_s*C_7eff;
 433     C_BR = -2*M_b*C_7eff;
 434     C_LLtot = C_9eff-m_WC.GetC10tilda()+C_LL;
 435     C_LRtot = C_9eff+m_WC.GetC10tilda()+C_LR;
 436     //std::cout << "Yld: " << Yld << "  C7eff: " << C_7eff << "  C_9eff: " << C_9eff << "  Diff7: " << C_7eff-m_WC.GetC7eff0() << "  Diff9: " << C_9eff-m_WC.GetC9tilda() << std::endl;
 437   }else if(m_HEPmodel=="-C7_SM"){
 438     C_LL=C_LR=C_RL=C_RR=C_LRLR=C_RLLR=C_LRRL=C_RLRL=C_T=C_TE=EvtComplex(0,0);
 439     Yld = m_WC.Yld(q2,k_psi,Gamma_psi,M_psi,2,m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),1./alpha);
 440     C_7eff = m_WC.GetC7eff0() + m_WC.C7b2sg(m_WC.GetStrongCouplingConst(),m_WC.GetEta(),m_WC.GetC2(),M_t,M_W) + m_omega*(m_WC.hzs(M_c/M_b,q2/M_b/M_b,M_b,M_b) + Yld);
 441     C_9eff = Yld + m_WC.C9efftilda(M_c/M_b,q2/M_b/M_b,m_WC.GetStrongCouplingConst(),m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),m_WC.GetC9tilda(),m_WC.GetRenormSchemePar());
 442     C_SL = +2*M_s*C_7eff;
 443     C_BR = +2*M_b*C_7eff;
 444     C_LLtot = C_9eff-m_WC.GetC10tilda()+C_LL;
 445     C_LRtot = C_9eff+m_WC.GetC10tilda()+C_LR;
 446     //std::cout << "Yld: " << Yld << "  C7eff: " << C_7eff << "  C_9eff: " << C_9eff << "  Diff7: " << C_7eff-m_WC.GetC7eff0() << "  Diff9: " << C_9eff-m_WC.GetC9tilda() << std::endl;
 447   }else if(m_HEPmodel=="SUSY-ChenGeng"){
 448     //report(WARNING,"EvtGen") << " WARNING: SUSY-ChenGeng model should be checked for correctness" << std::endl;
 449     C_LL=C_LR=C_RL=C_RR=C_LRLR=C_RLLR=C_LRRL=C_RLRL=C_T=C_TE=EvtComplex(0,0);
 450     EvtComplex d_u23LL = 0.1;
 451     EvtComplex d_u33RL = 0.65;
 452     EvtComplex d_d23LR = 0.03*exp(i1*EvtConst::pi*2/5);
 453     EvtComplex d_u23LR = -0.8*exp(i1*EvtConst::pi/4);
 454     EvtComplex C_7susy = -1.75*d_u23LL - 0.25*d_u23LR - 10.3*d_d23LR;
 455     EvtComplex C_9susy = 0.82*d_u23LR;
 456     EvtComplex C_10susy = -9.37*d_u23LR + 1.4*d_u23LR*d_u33RL + 2.7*d_u23LL;
 457     Yld = m_WC.Yld(q2,k_psi,Gamma_psi,M_psi,2,m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),1./alpha);
 458     C_7eff = m_WC.GetC7eff0() + C_7susy*pow(m_WC.GetEta(),16./23.) + m_WC.C7b2sg(m_WC.GetStrongCouplingConst(),m_WC.GetEta(),m_WC.GetC2(),M_t,M_W) + m_omega*(m_WC.hzs(M_c/M_b,q2/M_b/M_b,M_b,M_b) + Yld);
 459     C_9eff = Yld + m_WC.C9efftilda(M_c/M_b,q2/M_b/M_b,m_WC.GetStrongCouplingConst(),m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),m_WC.GetC9tilda()+C_9susy,m_WC.GetRenormSchemePar());
 460     C_SL = -2*M_s*C_7eff;
 461     C_BR = -2*M_b*C_7eff;
 462     C_LLtot = C_9eff-m_WC.GetC10tilda()-C_10susy+C_LL;
 463     C_LRtot = C_9eff+m_WC.GetC10tilda()+C_10susy+C_LR;
 464     //std::cout << "Yld: " << Yld << "  C7eff: " << C_7eff << "  C_9eff: " << C_9eff << "  Diff7: " << C_7eff-m_WC.GetC7eff0() << "  Diff9: " << C_9eff-m_WC.GetC9tilda() << std::endl;
 465   }else{
 466     report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll - unknown physics model: " << m_HEPmodel << " - this should never happen !" << std::endl;
 467     ::abort();
 468   }
 469
 470   A1 = (f_1T-g_1T)*C_SL/q2 + (f_1T+g_1T)*C_BR/q2 + 0.5*(f_1-g_1)*(C_LLtot+C_LRtot) + 0.5*(f_1+g_1)*(C_RL+C_RR);
 471   //std::cout << "f_1T: " << f_1T << "  g_1T: " << g_1T << "  C_SL: " << C_SL << "  C_BR: " << C_BR << "  f_1: " << f_1 << "  g_1: " << g_1 << "  C_LLtot: " << C_LLtot << "  C_LRtot: " << C_LRtot << "  C_RL: " << C_RL << "  C_RR: " << C_RR << std::endl;
 472   A2 = (f_2T-g_2T)*C_SL/q2 + (f_2T+g_2T)*C_BR/q2 + 0.5*(f_2-g_2)*(C_LLtot+C_LRtot) + 0.5*(f_2+g_2)*(C_RL+C_RR);
 473   A3 = (f_3T-g_3T)*C_SL/q2 + (f_3T+g_3T)*C_BR/q2 + 0.5*(f_3-g_3)*(C_LLtot+C_LRtot) + 0.5*(f_3+g_3)*(C_RL+C_RR);
 474
 475   B1 = (f_1T+g_1T)*C_SL/q2 + (f_1T-g_1T)*C_BR/q2 + 0.5*(f_1+g_1)*(C_LLtot+C_LRtot) + 0.5*(f_1-g_1)*(C_RL+C_RR);
 476   B2 = (f_2T+g_2T)*C_SL/q2 + (f_2T-g_2T)*C_BR/q2 + 0.5*(f_2+g_2)*(C_LLtot+C_LRtot) + 0.5*(f_2-g_2)*(C_RL+C_RR);
 477   B3 = (f_3T+g_3T)*C_SL/q2 + (f_3T-g_3T)*C_BR/q2 + 0.5*(f_3+g_3)*(C_LLtot+C_LRtot) + 0.5*(f_3-g_3)*(C_RL+C_RR);
 478
 479   D1 = 0.5*(C_RR-C_RL)*(f_1+g_1) + 0.5*(C_LRtot-C_LLtot)*(f_1-g_1);
 480   D2 = 0.5*(C_RR-C_RL)*(f_2+g_2) + 0.5*(C_LRtot-C_LLtot)*(f_2-g_2);
 481   D3 = 0.5*(C_RR-C_RL)*(f_3+g_3) + 0.5*(C_LRtot-C_LLtot)*(f_3-g_3);
 482
 483   E1 = 0.5*(C_RR-C_RL)*(f_1-g_1) + 0.5*(C_LRtot-C_LLtot)*(f_1+g_1);
 484   E2 = 0.5*(C_RR-C_RL)*(f_2-g_2) + 0.5*(C_LRtot-C_LLtot)*(f_2+g_2);
 485   E3 = 0.5*(C_RR-C_RL)*(f_3-g_3) + 0.5*(C_LRtot-C_LLtot)*(f_3+g_3);
 486
 487   N1 = (f_1*(M_Lb-M_L)+f_3*q2)/M_b*(C_LRLR+C_RLLR+C_LRRL+C_RLRL);  // Should be mLb - mL
 488   N2 = (f_1*(M_Lb-M_L)+f_3*q2)/M_b*(C_LRLR+C_RLLR-C_LRRL-C_RLRL);
 489
 490   H1 = (g_1*(M_Lb+M_L)-g_3*q2)/M_b*(C_LRLR-C_RLLR+C_LRRL-C_RLRL);
 491   H2 = (g_1*(M_Lb+M_L)-g_3*q2)/M_b*(C_LRLR-C_RLLR-C_LRRL+C_RLRL);
 492
 493
 494   for(i=0;i<4;i++){
 495     lbar_Gmu_l   [i/2][i%2] = EvtLeptonVCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2));
 496     lbar_GmuG5_l [i/2][i%2] = EvtLeptonACurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2));
 497     lbar_l       [i/2][i%2] = EvtLeptonSCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2));
 498     lbar_G5_l    [i/2][i%2] = EvtLeptonPCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2));
 499     lbar_Smunu_l [i/2][i%2] = EvtLeptonTCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2));
 500     lbar_ESmunu_l[i/2][i%2] = dual(EvtLeptonTCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2)));
 501   }
 502
 503   // TODO: polarization not yet introduced
 504   if(m_polarizationIntroduction=="SpinDensityMatrix"){
 505     //parent->setSpinDensityForward();
 506     parent__spParent[0]=parent->sp(0);
 507     parent__spParent[1]=parent->sp(1);
 508   }else if(m_polarizationIntroduction=="ModifiedSpinors"){
 509     parent__spParent[0]=parent->sp(0);
 510     parent__spParent[1]=parent->sp(1);
 511   }else if(m_polarizationIntroduction=="Unpolarized"){
 512     parent__spParent[0]=parent->sp(0);
 513     parent__spParent[1]=parent->sp(1);
 514   }else{
 515     report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll - unknown polarization: " << m_polarizationIntroduction << " - this should never happen !" << std::endl;
 516     ::abort();
 517   }
 518
 519   for(i=0;i<4;i++){
 520     hbar_GmuPlusG5_h   [i/2][i%2] = EvtLeptonVCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) + EvtLeptonACurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 521     hbar_GmuMinusG5_h  [i/2][i%2] = EvtLeptonVACurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 522     hbar_SmunuPlusG5_h [i/2][i%2] = EvtLeptonTCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) + EvtLeptonTG5Current(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 523     hbar_SmunuMinusG5_h[i/2][i%2] = EvtLeptonTCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) - EvtLeptonTG5Current(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 524     hbar_1PlusG5_h     [i/2][i%2] = EvtLeptonSCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) + EvtLeptonPCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 525     hbar_1MinusG5_h    [i/2][i%2] = EvtLeptonSCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) - EvtLeptonPCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 526     hbar_G5_h          [i/2][i%2] = EvtLeptonPCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 527     hbar_h             [i/2][i%2] = EvtLeptonSCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 528     hbar_Smunu_h       [i/2][i%2] = EvtLeptonTCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 529     hbar_Gmu_h         [i/2][i%2] = EvtLeptonVCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]);
 530   }
 531
 532   for(i=0;i<4;i++) for(j=0;j<4;j++) {
 533     //std::cout << "--------------------------------------------------" << std::endl;
 534     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 535     Matrix[j/2][j%2][i/2][i%2] += lbar_Gmu_l[i/2][i%2]    * (A1*hbar_GmuPlusG5_h[j/2][j%2]+B1*hbar_GmuMinusG5_h[j/2][j%2]);
 536     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 537     //std::cout << "A1: " << A1 << " B1: " << B1 << " lbar_Gmu_l: " << lbar_Gmu_l[i/2][i%2] <<
 538     //             " hbar_GmuPlusG5_h: " << hbar_GmuPlusG5_h[j/2][j%2] << " hbar_GmuMinusG5_h: " << hbar_GmuMinusG5_h[j/2][j%2] <<
 539     //             " sp1: " << parent->getDaug(1)->spParent(i/2) << " sp2: " << parent->getDaug(1)->spParent(i%2) << std::endl;
 540     Matrix[j/2][j%2][i/2][i%2] += lbar_Gmu_l[i/2][i%2]    * (i1*A2*(hbar_SmunuPlusG5_h[j/2][j%2].cont2(q_mu))+B2*(hbar_SmunuMinusG5_h[j/2][j%2].cont2(q_mu)));
 541     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 542     Matrix[j/2][j%2][i/2][i%2] += lbar_Gmu_l[i/2][i%2]    * ((A3*hbar_1PlusG5_h[j/2][j%2]+B3*hbar_1MinusG5_h[j/2][j%2])*q_mu);
 543     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 544     Matrix[j/2][j%2][i/2][i%2] += lbar_GmuG5_l[i/2][i%2]  * (D1*hbar_GmuPlusG5_h[j/2][j%2]+E1*hbar_GmuMinusG5_h[j/2][j%2]);
 545     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 546     Matrix[j/2][j%2][i/2][i%2] += lbar_GmuG5_l[i/2][i%2]  * (i1*D2*(hbar_SmunuPlusG5_h[j/2][j%2].cont2(q_mu))+E2*(hbar_SmunuMinusG5_h[j/2][j%2].cont2(q_mu)));
 547     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 548     Matrix[j/2][j%2][i/2][i%2] += lbar_GmuG5_l[i/2][i%2]  * ((D3*hbar_1PlusG5_h[j/2][j%2]+E3*hbar_1MinusG5_h[j/2][j%2])*q_mu);
 549     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 550     Matrix[j/2][j%2][i/2][i%2] += lbar_l[i/2][i%2]        * (N1*hbar_h[j/2][j%2]+H1*hbar_G5_h[j/2][j%2]);
 551     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 552     Matrix[j/2][j%2][i/2][i%2] += lbar_G5_l[i/2][i%2]     * (N2*hbar_h[j/2][j%2]+H2*hbar_G5_h[j/2][j%2]);
 553     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 554     Matrix[j/2][j%2][i/2][i%2] += cont(lbar_Smunu_l[i/2][i%2]  , 4*C_T*f_T*hbar_Smunu_h[j/2][j%2]);
 555     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 556     Matrix[j/2][j%2][i/2][i%2] +=
 557       cont(lbar_Smunu_l[i/2][i%2]  ,
 558            -4*C_T*f_TV*i1*(EvtGenFunctions::directProd(q_mu,hbar_Gmu_h[j/2][j%2])-
 559                            EvtGenFunctions::directProd(hbar_Gmu_h[j/2][j%2],q_mu)));
 560     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 561     Matrix[j/2][j%2][i/2][i%2] += cont(lbar_Smunu_l[i/2][i%2]  , -4*C_T*f_TS*i1*(EvtGenFunctions::directProd(P_mu,q_mu)-EvtGenFunctions::directProd(q_mu,P_mu))*hbar_h[j/2][j%2]);
 562     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 563     Matrix[j/2][j%2][i/2][i%2] += cont(lbar_ESmunu_l[i/2][i%2] , 4*C_TE*f_T*i1*hbar_Smunu_h[j/2][j%2]);
 564     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 565     Matrix[j/2][j%2][i/2][i%2] += cont(lbar_ESmunu_l[i/2][i%2] , 4*C_TE*f_TV*(EvtGenFunctions::directProd(q_mu,hbar_Gmu_h[j/2][j%2])-EvtGenFunctions::directProd(hbar_Gmu_h[j/2][j%2],q_mu)));
 566     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 567     Matrix[j/2][j%2][i/2][i%2] += cont(lbar_ESmunu_l[i/2][i%2] , 4*C_TE*f_TS*(EvtGenFunctions::directProd(P_mu,q_mu)-EvtGenFunctions::directProd(q_mu,P_mu))*hbar_h[j/2][j%2]);
 568     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 569     //Matrix[j/2][j%2][i/2][i%2] *= G_F*alpha/4/sqrt(2)/EvtConst::pi*V_tb*conj(V_ts);
 570     //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
 571     //std::cout << "--------------------------------------------------" << std::endl;
 572     spins[0]=j%2;
 573     spins[1]=j/2;
 574     spins[2]=i/2;
 575     spins[3]=i%2;
 576     amp->vertex(spins,Matrix[j/2][j%2][i/2][i%2]);
 577   }
 578
 579   //std::cout << "==================================================" << std::endl;
 580   //std::cout << "Lambda_b0:  " << parent->getP4Restframe() << std::endl;
 581   //std::cout << "Lambda0:  " << parent->getDaug(0)->getP4() << std::endl;
 582   //std::cout << "mu-:  " << parent->getDaug(1)->getP4() << std::endl;
 583   //std::cout << "mu+:  " << parent->getDaug(2)->getP4() << std::endl;
 584   //std::cout << "P_mu:  " << P_mu << std::endl;
 585   //std::cout << "q_mu:  " << q_mu << std::endl;
 586   //std::cout << "q2:  " << q2 << std::endl;
 587   //std::cout << "==================================================" << std::endl;
 588
 589   return;
 590 }
 591
 592 EvtTensor4C EvtLb2Lll::EvtLeptonTG5Current(const EvtDiracSpinor &d,const EvtDiracSpinor &dp){
 593 // <u|sigma^munu*gamma^5|v>
 594
 595   EvtTensor4C temp;
 596   temp.zero();
 597   EvtComplex i2(0,0.5);
 598
 599   static EvtGammaMatrix mat01=EvtGammaMatrix::g0()*
 600     (EvtGammaMatrix::g0()*EvtGammaMatrix::g1()-
 601      EvtGammaMatrix::g1()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5();
 602   static EvtGammaMatrix mat02=EvtGammaMatrix::g0()*
 603     (EvtGammaMatrix::g0()*EvtGammaMatrix::g2()-
 604      EvtGammaMatrix::g2()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5();
 605   static EvtGammaMatrix mat03=EvtGammaMatrix::g0()*
 606     (EvtGammaMatrix::g0()*EvtGammaMatrix::g3()-
 607      EvtGammaMatrix::g3()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5();
 608   static EvtGammaMatrix mat12=EvtGammaMatrix::g0()*
 609     (EvtGammaMatrix::g1()*EvtGammaMatrix::g2()-
 610      EvtGammaMatrix::g2()*EvtGammaMatrix::g1())*EvtGammaMatrix::g5();
 611   static EvtGammaMatrix mat13=EvtGammaMatrix::g0()*
 612     (EvtGammaMatrix::g1()*EvtGammaMatrix::g3()-
 613      EvtGammaMatrix::g3()*EvtGammaMatrix::g1())*EvtGammaMatrix::g5();
 614   static EvtGammaMatrix mat23=EvtGammaMatrix::g0()*
 615     (EvtGammaMatrix::g2()*EvtGammaMatrix::g3()-
 616      EvtGammaMatrix::g3()*EvtGammaMatrix::g2())*EvtGammaMatrix::g5();
 617
 618   temp.set(0,1,i2*(d*(mat01*dp)));
 619   temp.set(1,0,-temp.get(0,1));
 620
 621   temp.set(0,2,i2*(d*(mat02*dp)));
 622   temp.set(2,0,-temp.get(0,2));
 623
 624   temp.set(0,3,i2*(d*(mat03*dp)));
 625   temp.set(3,0,-temp.get(0,3));
 626
 627   temp.set(1,2,i2*(d*(mat12*dp)));
 628   temp.set(2,1,-temp.get(1,2));
 629
 630   temp.set(1,3,i2*(d*(mat13*dp)));
 631   temp.set(3,1,-temp.get(1,3));
 632
 633   temp.set(2,3,i2*(d*(mat23*dp)));
 634   temp.set(3,2,-temp.get(2,3));
 635
 636   return temp;
 637 }