[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtbTosllAmp.cxx

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtbTosllAmp.cc
//
// Description: Routine to implement semileptonic decays to pseudo-scalar
//              mesons.
//
// Modification history:
//
//    DJL       April 17,1998       Module created
//
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenModels/EvtbTosllAmp.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtScalarParticle.hh"
#include "EvtGenBase/EvtVectorParticle.hh"
#include "EvtGenBase/EvtDiLog.hh"

double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson, 
				  EvtId lepton1, EvtId lepton2,
				  EvtbTosllFF *FormFactors,
				  double& poleSize) {

  //This routine takes the arguements parent, meson, and lepton
  //number, and a form factor model, and returns a maximum
  //probability for this semileptonic form factor model.  A
  //brute force method is used.  The 2D cos theta lepton and
  //q2 phase space is probed.

  //Start by declaring a particle at rest.

  //It only makes sense to have a scalar parent.  For now. 
  //This should be generalized later.

  EvtScalarParticle *scalar_part;
  EvtParticle *root_part;

  scalar_part=new EvtScalarParticle;

  //cludge to avoid generating random numbers!
  scalar_part->noLifeTime();

  EvtVector4R p_init;
  
  p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
  scalar_part->init(parent,p_init);
  root_part=(EvtParticle *)scalar_part;
  root_part->setDiagonalSpinDensity();      

  EvtParticle *daughter, *lep1, *lep2;
  
  EvtAmp amp;

  EvtId listdaug[3];
  listdaug[0] = meson;
  listdaug[1] = lepton1;
  listdaug[2] = lepton2;

  amp.init(parent,3,listdaug);

  root_part->makeDaughters(3,listdaug);
  daughter=root_part->getDaug(0);
  lep1=root_part->getDaug(1);
  lep2=root_part->getDaug(2);

  //cludge to avoid generating random numbers!
  daughter->noLifeTime();
  lep1->noLifeTime();
  lep2->noLifeTime();


  //Initial particle is unpolarized, well it is a scalar so it is 
  //trivial
  EvtSpinDensity rho;
  rho.setDiag(root_part->getSpinStates());
  
  double mass[3];
  
  double m = root_part->mass();
  
  EvtVector4R p4meson, p4lepton1, p4lepton2, p4w;
  double q2max;

  double q2, elepton, plepton;
  int i,j;
  double erho,prho,costl;

  double maxfoundprob = 0.0;
  double prob = -10.0;
  int massiter;

  double maxpole=0;

  for (massiter=0;massiter<3;massiter++){

    mass[0] = EvtPDL::getMeanMass(meson);
    mass[1] = EvtPDL::getMeanMass(lepton1);
    mass[2] = EvtPDL::getMeanMass(lepton2);
    if ( massiter==1 ) {
      mass[0] = EvtPDL::getMinMass(meson);
    }
    if ( massiter==2 ) {
      mass[0] = EvtPDL::getMaxMass(meson);
      if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001; 
    }

    q2max = (m-mass[0])*(m-mass[0]);
    
    //loop over q2
    //cout << "m " << m << "mass[0] " << mass[0] << " q2max "<< q2max << endl;
    for (i=0;i<25;i++) {
      //want to avoid picking up the tail of the photon propagator
      q2 = ((i+1.5)*q2max)/26.0;

      if (i==0) q2=4*(mass[1]*mass[1]);

      erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
      
      prho = sqrt(erho*erho-mass[0]*mass[0]);
      
      p4meson.set(erho,0.0,0.0,-1.0*prho);
      p4w.set(m-erho,0.0,0.0,prho);
      
      //This is in the W rest frame
      elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
      plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
      
      double probctl[3];

      for (j=0;j<3;j++) {
	
        costl = 0.99*(j - 1.0);
	
	//These are in the W rest frame. Need to boost out into
	//the B frame.
        p4lepton1.set(elepton,0.0,
		  plepton*sqrt(1.0-costl*costl),plepton*costl);
        p4lepton2.set(elepton,0.0,
		 -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);

	EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
        p4lepton1=boostTo(p4lepton1,boost);
        p4lepton2=boostTo(p4lepton2,boost);

	//Now initialize the daughters...

        daughter->init(meson,p4meson);
        lep1->init(lepton1,p4lepton1);
        lep2->init(lepton2,p4lepton2);

        CalcAmp(root_part,amp,FormFactors);

	//Now find the probability at this q2 and cos theta lepton point
        //and compare to maxfoundprob.

	//Do a little magic to get the probability!!

	//cout <<"amp:"<<amp.getSpinDensity()<<endl;

	prob = rho.normalizedProb(amp.getSpinDensity());

	//cout << "prob:"<<q2<<" "<<costl<<" "<<prob<<endl;

	probctl[j]=prob;
      }

      //probclt contains prob at ctl=-1,0,1.
      //prob=a+b*ctl+c*ctl^2

      double a=probctl[1];
      double b=0.5*(probctl[2]-probctl[0]);
      double c=0.5*(probctl[2]+probctl[0])-probctl[1];

      prob=probctl[0];
      if (probctl[1]>prob) prob=probctl[1];
      if (probctl[2]>prob) prob=probctl[2];

      if (fabs(c)>1e-20){
	double ctlx=-0.5*b/c;
	if (fabs(ctlx)<1.0){
	  double probtmp=a+b*ctlx+c*ctlx*ctlx;
	  if (probtmp>prob) prob=probtmp;
	} 

      }

      //report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
      //			    << probctl[0]<<" "
      //			    << probctl[1]<<" "
      //			    << probctl[2]<<endl;

      if (i==0) {
	maxpole=prob;
	continue;
      }

      if ( prob > maxfoundprob ) {
	maxfoundprob = prob; 
      }

      //cout << "q2,maxfoundprob:"<<q2<<" "<<maxfoundprob<<endl;

    }
    if ( EvtPDL::getWidth(meson) <= 0.0 ) {
      //if the particle is narrow dont bother with changing the mass.
      massiter = 4;
    }

  }

  root_part->deleteTree();  

  poleSize=0.04*(maxpole/maxfoundprob)*4*(mass[1]*mass[1]);

  //poleSize=0.002;

  //cout <<"maxfoundprob,maxpole,poleSize:"<<maxfoundprob<<" "
  //     <<maxpole<<" "<<poleSize<<endl;

  maxfoundprob *=1.15;

  return maxfoundprob;
  
}


EvtComplex EvtbTosllAmp::GetC7Eff(double q2, bool nnlo) 
{

  if (!nnlo) return -0.313;
  double mbeff = 4.8;
  double shat = q2/mbeff/mbeff;
  double logshat;
  logshat = log(shat);
  
  double muscale;
  muscale = 2.5;
  double alphas;
  alphas = 0.267;
  double A7;
  A7 = -0.353 + 0.023;
  double A8;
  A8 = -0.164;
  double A9;
  A9 = 4.287 + (-0.218);
  double A10;
  A10 = -4.592 + 0.379;
  double C1;
  C1 = -0.697;
  double C2;
  C2 = 1.046;
  double T9;
  T9 = 0.114 + 0.280;
  double U9;
  U9 = 0.045 + 0.023;
  double W9;
  W9 = 0.044 + 0.016;
  
  double Lmu;
  Lmu = log(muscale/mbeff);

  EvtComplex uniti(0.0,1.0);

  EvtComplex c7eff;
  if (shat > 0.25)
  { 
   c7eff = A7;
   return c7eff;
  }


  // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
  muscale = 5.0;
  alphas = 0.215;
  A7 = -0.312 + 0.008;
  A8 = -0.148;
  A9 = 4.174 + (-0.035);
  A10 = -4.592 + 0.379;
  C1 = -0.487;
  C2 = 1.024;
  T9 = 0.374 + 0.252;
  U9 = 0.033 + 0.015;
  W9 = 0.032 + 0.012;
  Lmu = log(muscale/mbeff);

  EvtComplex F71;
  EvtComplex f71;
  EvtComplex k7100(-0.68192,-0.074998);
  EvtComplex k7101(0.0,0.0);
  EvtComplex k7110(-0.23935,-0.12289);
  EvtComplex k7111(0.0027424,0.019676);
  EvtComplex k7120(-0.0018555,-0.175);
  EvtComplex k7121(0.022864,0.011456);
  EvtComplex k7130(0.28248,-0.12783);
  EvtComplex k7131(0.029027,-0.0082265);
  f71 = k7100 + k7101*logshat + shat*(k7110 + k7111*logshat) +
        shat*shat*(k7120 + k7121*logshat) + 
        shat*shat*shat*(k7130 + k7131*logshat); 
  F71 = (-208.0/243.0)*Lmu + f71;

  EvtComplex F72;
  EvtComplex f72;
  EvtComplex k7200(4.0915,0.44999);
  EvtComplex k7201(0.0,0.0);
  EvtComplex k7210(1.4361,0.73732);
  EvtComplex k7211(-0.016454,-0.11806);
  EvtComplex k7220(0.011133,1.05);
  EvtComplex k7221(-0.13718,-0.068733);
  EvtComplex k7230(-1.6949,0.76698);
  EvtComplex k7231(-0.17416,0.049359);
  f72 = k7200 + k7201*logshat + shat*(k7210 + k7211*logshat) +
        shat*shat*(k7220 + k7221*logshat) + 
        shat*shat*shat*(k7230 + k7231*logshat); 
  F72 = (416.0/81.0)*Lmu + f72;
  
  EvtComplex F78;
  F78 = (-32.0/9.0)*Lmu + 8.0*EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0) 
        + (-8.0*EvtConst::pi/9.0)*uniti +
        (4.0/3.0*EvtConst::pi*EvtConst::pi - 40.0/3.0)*shat +
        (32.0*EvtConst::pi*EvtConst::pi/9.0 - 316.0/9.0)*shat*shat +
        (200.0*EvtConst::pi*EvtConst::pi/27.0 - 658.0/9.0)*shat*shat*shat +
    (-8.0*logshat/9.0)*(shat + shat*shat + shat*shat*shat);
        
  c7eff = A7 - alphas/(4.0*EvtConst::pi)*(C1*F71 + C2*F72 + A8*F78);

  return c7eff;
}


EvtComplex EvtbTosllAmp::GetC9Eff(double q2, bool nnlo, bool btod) 
{

  if (!nnlo) return 4.344;
  double mbeff = 4.8;
  double shat = q2/mbeff/mbeff;
  double logshat;
  logshat = log(shat);
  double mchat = 0.29;

  
  double muscale;
  muscale = 2.5;
  double alphas;
  alphas = 0.267;
  double A7;
  A7 = -0.353 + 0.023;
  double A8;
  A8 = -0.164;
  double A9;
  A9 = 4.287 + (-0.218);
  double A10;
  A10 = -4.592 + 0.379;
  double C1;
  C1 = -0.697;
  double C2;
  C2 = 1.046;
  double T9;
  T9 = 0.114 + 0.280;
  double U9;
  U9 = 0.045 + 0.023;
  double W9;
  W9 = 0.044 + 0.016;
  
  double Lmu;
  Lmu = log(muscale/mbeff);


  EvtComplex uniti(0.0,1.0);

  EvtComplex hc;
  double xarg;
  xarg = 4.0*mchat/shat;
  hc = -4.0/9.0*log(mchat*mchat) + 8.0/27.0 + 4.0*xarg/9.0;

if (xarg < 1.0)
  {
    hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
      (log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
       uniti*EvtConst::pi);
  }
  else
  {
    hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
      2.0*atan(1.0/sqrt(xarg - 1.0));
  }
                                                                                                                                                             
  EvtComplex h1;
  xarg = 4.0/shat;
  h1 = 8.0/27.0 + 4.0*xarg/9.0;
  if (xarg < 1.0)
  {
    h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
      (log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
       uniti*EvtConst::pi);
  }
  else
  {
    h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
      2.0*atan(1.0/sqrt(xarg - 1.0));
  }


  EvtComplex h0;
  h0 = 8.0/27.0 - 4.0*log(2.0)/9.0 + 4.0*uniti*EvtConst::pi/9.0;


  // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
  // h(\hat m_u^2, hat s))
  EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
  EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
  EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
  EvtComplex Vtb(1.0,0.0);

  EvtComplex Xd;
  Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);


  EvtComplex c9eff=4.344;
  if (shat > 0.25)
  { 
   c9eff =  A9 + T9*hc + U9*h1 + W9*h0;
   if (btod)
   {
    c9eff += Xd; 
   }

   return c9eff;
  }

  // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
  muscale = 5.0;
  alphas = 0.215;
  A9 = 4.174 + (-0.035);
  C1 = -0.487;
  C2 = 1.024;
  A8 = -0.148;
  T9 = 0.374 + 0.252;
  U9 = 0.033 + 0.015;
  W9 = 0.032 + 0.012;
  Lmu = log(muscale/mbeff);

  EvtComplex F91;
  EvtComplex f91;
  EvtComplex k9100(-11.973,0.16371);
  EvtComplex k9101(-0.081271,-0.059691);
  EvtComplex k9110(-28.432,-0.25044);
  EvtComplex k9111(-0.040243,0.016442);
  EvtComplex k9120(-57.114,-0.86486);
  EvtComplex k9121(-0.035191,0.027909);
  EvtComplex k9130(-128.8,-2.5243);
  EvtComplex k9131(-0.017587,0.050639);
  f91 = k9100 + k9101*logshat + shat*(k9110 + k9111*logshat) +
        shat*shat*(k9120 + k9121*logshat) + 
        shat*shat*shat*(k9130 + k9131*logshat); 
  F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0 
         + 64.0/27.0*log(mchat))*Lmu - 16.0*Lmu*logshat/243.0 +
        (16.0/1215.0 - 32.0/135.0/mchat/mchat)*Lmu*shat +
        (4.0/2835.0 - 8.0/315.0/mchat/mchat/mchat/mchat)*Lmu*shat*shat +
    (16.0/76545.0 - 32.0/8505.0/mchat/mchat/mchat/mchat/mchat/mchat)*
    Lmu*shat*shat*shat -256.0*Lmu*Lmu/243.0 + f91;

  EvtComplex F92;
  EvtComplex f92;
  EvtComplex k9200(6.6338,-0.98225);
  EvtComplex k9201(0.48763,0.35815);
  EvtComplex k9210(3.3585,1.5026);
  EvtComplex k9211(0.24146,-0.098649);
  EvtComplex k9220(-1.1906,5.1892);
  EvtComplex k9221(0.21115,-0.16745);
  EvtComplex k9230(-17.12,15.146);
  EvtComplex k9231(0.10552,-0.30383);
  f92 = k9200 + k9201*logshat + shat*(k9210 + k9211*logshat) +
        shat*shat*(k9220 + k9221*logshat) + 
        shat*shat*shat*(k9230 + k9231*logshat); 
  F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0 
         - 128.0/9.0*log(mchat))*Lmu + 32.0*Lmu*logshat/81.0 +
        (-32.0/405.0 + 64.0/45.0/mchat/mchat)*Lmu*shat +
        (-8.0/945.0 + 16.0/105.0/mchat/mchat/mchat/mchat)*Lmu*shat*shat +
    (-32.0/25515.0 + 64.0/2835.0/mchat/mchat/mchat/mchat/mchat/mchat)*
    Lmu*shat*shat*shat + 512.0*Lmu*Lmu/81.0 + f92;
  
  EvtComplex F98;
  F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 + 
        (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*shat +
        (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*shat*shat +
    (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*shat*shat*shat +
        16.0*logshat/9.0*(1.0 + shat + shat*shat + shat*shat*shat);

  Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);

  c9eff = A9 + T9*hc + U9*h1 + W9*h0 -             
    alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98);
  if (btod)
  {
   c9eff += Xd; 
  }

  return c9eff;
}

EvtComplex EvtbTosllAmp::GetC10Eff(double /*q2*/, bool nnlo) 
{

  if (!nnlo) return -4.669;
  double A10;
  A10 = -4.592 + 0.379;

  EvtComplex c10eff;
  c10eff = A10;

  return c10eff;
}

double EvtbTosllAmp::dGdsProb(double mb, double ms, double ml,
                                double s)
{
  // Compute the decay probability density function given a value of s
  // according to Ali's paper


  double delta, lambda, prob;
  double f1, f2, f3, f4;
  double msh, mlh, sh;

  mlh = ml / mb;
  msh = ms / mb;
  sh  = s  / (mb*mb);

  EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
  EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
  EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);

  double alphas = 0.119/
     (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
  double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
                 - 2.0/3.0*log(sh)*log(1.0-sh)
                 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
                 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
                 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
                 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
  double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
  double omega7 = -8.0/3.0*log(4.8/mb)
                  -4.0/3.0*EvtDiLog::DiLog(sh) 
                  -2.0/9.0*EvtConst::pi*EvtConst::pi
                  -2.0/3.0*log(sh)*log(1.0-sh)
                  -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0 
    -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
    -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
  double eta7 = 1.0 + alphas*omega7/EvtConst::pi;

  double omega79 = -4.0/3.0*log(4.8/mb)
                   -4.0/3.0*EvtDiLog::DiLog(sh) 
                   -2.0/9.0*EvtConst::pi*EvtConst::pi
                   -2.0/3.0*log(sh)*log(1.0-sh)
                   -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
                   -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2) 
                   +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
  double eta79 = 1.0 + alphas*omega79/EvtConst::pi;

  double c7c9 = abs(c7eff)*real(c9eff);
  c7c9 *= pow(eta79,2); 
  double c7c7 = pow(abs(c7eff),2);
  c7c7 *= pow(eta7,2); 

  double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
  c9c9plusc10c10 *= pow(eta9,2);
  double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
  c9c9minusc10c10 *= pow(eta9,2);
 
  lambda = 1.0 + sh*sh + pow(msh,4) - 2.0*(sh + sh*msh*msh + msh*msh);

  f1 = pow(1.0-msh*msh,2) - sh*(1.0 + msh*msh);
  f2 = 2.0*(1.0 + msh*msh) * pow(1.0-msh*msh,2)
       - sh*(1.0 + 14.0*msh*msh + pow(msh,4)) - sh*sh*(1.0 + msh*msh);
  f3 = pow(1.0-msh*msh,2) + sh*(1.0 + msh*msh) - 2.0*sh*sh
       + lambda*2.0*mlh*mlh/sh;
  f4 = 1.0 - sh + msh*msh;

  delta = (  12.0*c7c9*f1 + 4.0*c7c7*f2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
            + c9c9plusc10c10*f3 
            + 6.0*mlh*mlh*c9c9minusc10c10*f4;

  prob =  sqrt(lambda*(1.0 - 4.0*mlh*mlh/sh)) * delta;

  return prob;
}

double EvtbTosllAmp::dGdsdupProb(double mb, double ms, double ml,
                                   double s,  double u)
{
  // Compute the decay probability density function given a value of s and u
  // according to Ali's paper

  double prob;
  double f1sp, f2sp, f3sp;

  double sh = s / (mb*mb);

  EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
  EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
  EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);

  double alphas = 0.119/
     (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
  double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
                 - 2.0/3.0*log(sh)*log(1.0-sh)
                 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
                 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
                 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
                 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
  double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
  double omega7 = -8.0/3.0*log(4.8/mb)
                  -4.0/3.0*EvtDiLog::DiLog(sh) 
                  -2.0/9.0*EvtConst::pi*EvtConst::pi
                  -2.0/3.0*log(sh)*log(1.0-sh)
                  -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0 
    -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
    -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
  double eta7 = 1.0 + alphas*omega7/EvtConst::pi;

  double omega79 = -4.0/3.0*log(4.8/mb)
                   -4.0/3.0*EvtDiLog::DiLog(sh) 
                   -2.0/9.0*EvtConst::pi*EvtConst::pi
                   -2.0/3.0*log(sh)*log(1.0-sh)
                   -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
                   -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2) 
                   +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
  double eta79 = 1.0 + alphas*omega79/EvtConst::pi;

  double c7c9 = abs(c7eff)*real(c9eff);
  c7c9 *= pow(eta79,2); 
  double c7c7 = pow(abs(c7eff),2);
  c7c7 *= pow(eta7,2); 

  double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
  c9c9plusc10c10 *= pow(eta9,2);
  double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
  c9c9minusc10c10 *= pow(eta9,2);
  double c7c10 = abs(c7eff)*real(c10eff);
  c7c10 *= eta7; c7c10 *= eta9;
  double c9c10 = real(c9eff)*real(c10eff);
  c9c10 *= pow(eta9,2); 

  f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10 
         + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb))
         - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s
    // kludged mass term
         *(1.0 + 2.0*ml*ml/s)
         - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9
    // kludged mass term
         *(1.0 + 2.0*ml*ml/s);

  f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10;
  f3sp = - (c9c9plusc10c10)
         + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s
    // kludged mass term
         *(1.0 + 2.0*ml*ml/s);

  prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3);

  return prob;
}
Commit	Line	Data
da0e9ce3	1	//--------------------------------------------------------------------------
	2	//
	3	// Environment:
	4	// This software is part of the EvtGen package developed jointly
	5	// for the BaBar and CLEO collaborations. If you use all or part
	6	// of it, please give an appropriate acknowledgement.
	7	//
	8	// Copyright Information: See EvtGen/COPYRIGHT
	9	// Copyright (C) 1998 Caltech, UCSB
	10	//
	11	// Module: EvtbTosllAmp.cc
	12	//
	13	// Description: Routine to implement semileptonic decays to pseudo-scalar
	14	// mesons.
	15	//
	16	// Modification history:
	17	//
	18	// DJL April 17,1998 Module created
	19	//
	20	//------------------------------------------------------------------------
	21	//
	22	#include "EvtGenBase/EvtPatches.hh"
	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/EvtVector4C.hh"
	29	#include "EvtGenBase/EvtTensor4C.hh"
	30	#include "EvtGenBase/EvtDiracSpinor.hh"
	31	#include "EvtGenModels/EvtbTosllAmp.hh"
	32	#include "EvtGenBase/EvtId.hh"
	33	#include "EvtGenBase/EvtAmp.hh"
	34	#include "EvtGenBase/EvtScalarParticle.hh"
	35	#include "EvtGenBase/EvtVectorParticle.hh"
	36	#include "EvtGenBase/EvtDiLog.hh"
	37
	38	double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson,
	39	EvtId lepton1, EvtId lepton2,
	40	EvtbTosllFF *FormFactors,
	41	double& poleSize) {
	42
	43	//This routine takes the arguements parent, meson, and lepton
	44	//number, and a form factor model, and returns a maximum
	45	//probability for this semileptonic form factor model. A
	46	//brute force method is used. The 2D cos theta lepton and
	47	//q2 phase space is probed.
	48
	49	//Start by declaring a particle at rest.
	50
	51	//It only makes sense to have a scalar parent. For now.
	52	//This should be generalized later.
	53
	54	EvtScalarParticle *scalar_part;
	55	EvtParticle *root_part;
	56
	57	scalar_part=new EvtScalarParticle;
	58
	59	//cludge to avoid generating random numbers!
	60	scalar_part->noLifeTime();
	61
	62	EvtVector4R p_init;
	63
	64	p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
65	scalar_part->init(parent,p_init);
66	root_part=(EvtParticle *)scalar_part;
67	root_part->setDiagonalSpinDensity();
68
69	EvtParticle daughter, lep1, *lep2;
70
71	EvtAmp amp;
72
73	EvtId listdaug[3];
74	listdaug[0] = meson;
75	listdaug[1] = lepton1;
76	listdaug[2] = lepton2;
77
78	amp.init(parent,3,listdaug);
79
80	root_part->makeDaughters(3,listdaug);
81	daughter=root_part->getDaug(0);
82	lep1=root_part->getDaug(1);
83	lep2=root_part->getDaug(2);
84
85	//cludge to avoid generating random numbers!
86	daughter->noLifeTime();
87	lep1->noLifeTime();
88	lep2->noLifeTime();
89
90
91	//Initial particle is unpolarized, well it is a scalar so it is
92	//trivial
93	EvtSpinDensity rho;
94	rho.setDiag(root_part->getSpinStates());
95
96	double mass[3];
97
98	double m = root_part->mass();
99
100	EvtVector4R p4meson, p4lepton1, p4lepton2, p4w;
101	double q2max;
102
103	double q2, elepton, plepton;
104	int i,j;
105	double erho,prho,costl;
106
107	double maxfoundprob = 0.0;
108	double prob = -10.0;
109	int massiter;
110
111	double maxpole=0;
112
113	for (massiter=0;massiter<3;massiter++){
114
115	mass[0] = EvtPDL::getMeanMass(meson);
116	mass[1] = EvtPDL::getMeanMass(lepton1);
117	mass[2] = EvtPDL::getMeanMass(lepton2);
118	if ( massiter==1 ) {
119	mass[0] = EvtPDL::getMinMass(meson);
120	}
121	if ( massiter==2 ) {
122	mass[0] = EvtPDL::getMaxMass(meson);
123	if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
124	}
125
126	q2max = (m-mass[0])*(m-mass[0]);
127
128	//loop over q2
129	//cout << "m " << m << "mass[0] " << mass[0] << " q2max "<< q2max << endl;
130	for (i=0;i<25;i++) {
131	//want to avoid picking up the tail of the photon propagator
132	q2 = ((i+1.5)*q2max)/26.0;
133
134	if (i==0) q2=4(mass[1]mass[1]);
135
136	erho = ( mm + mass[0]mass[0] - q2 )/(2.0*m);
137
138	prho = sqrt(erhoerho-mass[0]mass[0]);
139
140	p4meson.set(erho,0.0,0.0,-1.0*prho);
141	p4w.set(m-erho,0.0,0.0,prho);
142
143	//This is in the W rest frame
144	elepton = (q2+mass[1]mass[1])/(2.0sqrt(q2));
145	plepton = sqrt(eleptonelepton-mass[1]mass[1]);
146
147	double probctl[3];
148
149	for (j=0;j<3;j++) {
150
151	costl = 0.99*(j - 1.0);
152
153	//These are in the W rest frame. Need to boost out into
154	//the B frame.
155	p4lepton1.set(elepton,0.0,
156	pleptonsqrt(1.0-costlcostl),plepton*costl);
157	p4lepton2.set(elepton,0.0,
158	-1.0pleptonsqrt(1.0-costlcostl),-1.0plepton*costl);
159
160	EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
161	p4lepton1=boostTo(p4lepton1,boost);
162	p4lepton2=boostTo(p4lepton2,boost);
163
164	//Now initialize the daughters...
165
166	daughter->init(meson,p4meson);
167	lep1->init(lepton1,p4lepton1);
168	lep2->init(lepton2,p4lepton2);
169
170	CalcAmp(root_part,amp,FormFactors);
171
172	//Now find the probability at this q2 and cos theta lepton point
173	//and compare to maxfoundprob.
174
175	//Do a little magic to get the probability!!
176
177	//cout <<"amp:"<<amp.getSpinDensity()<<endl;
178
179	prob = rho.normalizedProb(amp.getSpinDensity());
180
181	//cout << "prob:"<<q2<<" "<<costl<<" "<<prob<<endl;
182
183	probctl[j]=prob;
184	}
185
186	//probclt contains prob at ctl=-1,0,1.
187	//prob=a+bctl+cctl^2
188
189	double a=probctl[1];
190	double b=0.5*(probctl[2]-probctl[0]);
191	double c=0.5*(probctl[2]+probctl[0])-probctl[1];
192
193	prob=probctl[0];
194	if (probctl[1]>prob) prob=probctl[1];
195	if (probctl[2]>prob) prob=probctl[2];
196
197	if (fabs(c)>1e-20){
198	double ctlx=-0.5*b/c;
199	if (fabs(ctlx)<1.0){
200	double probtmp=a+bctlx+cctlx*ctlx;
201	if (probtmp>prob) prob=probtmp;
202	}
203
204	}
205
206	//report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
207	// << probctl[0]<<" "
208	// << probctl[1]<<" "
209	// << probctl[2]<<endl;
210
211	if (i==0) {
212	maxpole=prob;
213	continue;
214	}
215
216	if ( prob > maxfoundprob ) {
217	maxfoundprob = prob;
218	}
219
220	//cout << "q2,maxfoundprob:"<<q2<<" "<<maxfoundprob<<endl;
221
222	}
223	if ( EvtPDL::getWidth(meson) <= 0.0 ) {
224	//if the particle is narrow dont bother with changing the mass.
225	massiter = 4;
226	}
227
228	}
229
230	root_part->deleteTree();
231
232	poleSize=0.04(maxpole/maxfoundprob)4(mass[1]mass[1]);
233
234	//poleSize=0.002;
235
236	//cout <<"maxfoundprob,maxpole,poleSize:"<<maxfoundprob<<" "
237	// <<maxpole<<" "<<poleSize<<endl;
238
239	maxfoundprob *=1.15;
240
241	return maxfoundprob;
242
243	}
244
245
246	EvtComplex EvtbTosllAmp::GetC7Eff(double q2, bool nnlo)
247	{
248
249	if (!nnlo) return -0.313;
250	double mbeff = 4.8;
251	double shat = q2/mbeff/mbeff;
252	double logshat;
253	logshat = log(shat);
254
255	double muscale;
256	muscale = 2.5;
257	double alphas;
258	alphas = 0.267;
259	double A7;
260	A7 = -0.353 + 0.023;
261	double A8;
262	A8 = -0.164;
263	double A9;
264	A9 = 4.287 + (-0.218);
265	double A10;
266	A10 = -4.592 + 0.379;
267	double C1;
268	C1 = -0.697;
269	double C2;
270	C2 = 1.046;
271	double T9;
272	T9 = 0.114 + 0.280;
273	double U9;
274	U9 = 0.045 + 0.023;
275	double W9;
276	W9 = 0.044 + 0.016;
277
278	double Lmu;
279	Lmu = log(muscale/mbeff);
280
281	EvtComplex uniti(0.0,1.0);
282
283	EvtComplex c7eff;
284	if (shat > 0.25)
285	{
286	c7eff = A7;
287	return c7eff;
288	}
289
290
291
292
293	// change energy scale to 5.0 for full NNLO calculation below shat = 0.25
294	muscale = 5.0;
295	alphas = 0.215;
296	A7 = -0.312 + 0.008;
297	A8 = -0.148;
298	A9 = 4.174 + (-0.035);
299	A10 = -4.592 + 0.379;
300	C1 = -0.487;
301	C2 = 1.024;
302	T9 = 0.374 + 0.252;
303	U9 = 0.033 + 0.015;
304	W9 = 0.032 + 0.012;
305	Lmu = log(muscale/mbeff);
306
307	EvtComplex F71;
308	EvtComplex f71;
309	EvtComplex k7100(-0.68192,-0.074998);
310	EvtComplex k7101(0.0,0.0);
311	EvtComplex k7110(-0.23935,-0.12289);
312	EvtComplex k7111(0.0027424,0.019676);
313	EvtComplex k7120(-0.0018555,-0.175);
314	EvtComplex k7121(0.022864,0.011456);
315	EvtComplex k7130(0.28248,-0.12783);
316	EvtComplex k7131(0.029027,-0.0082265);
317	f71 = k7100 + k7101logshat + shat(k7110 + k7111*logshat) +
318	shatshat(k7120 + k7121*logshat) +
319	shatshatshat(k7130 + k7131logshat);
320	F71 = (-208.0/243.0)*Lmu + f71;
321
322	EvtComplex F72;
323	EvtComplex f72;
324	EvtComplex k7200(4.0915,0.44999);
325	EvtComplex k7201(0.0,0.0);
326	EvtComplex k7210(1.4361,0.73732);
327	EvtComplex k7211(-0.016454,-0.11806);
328	EvtComplex k7220(0.011133,1.05);
329	EvtComplex k7221(-0.13718,-0.068733);
330	EvtComplex k7230(-1.6949,0.76698);
331	EvtComplex k7231(-0.17416,0.049359);
332	f72 = k7200 + k7201logshat + shat(k7210 + k7211*logshat) +
333	shatshat(k7220 + k7221*logshat) +
334	shatshatshat(k7230 + k7231logshat);
335	F72 = (416.0/81.0)*Lmu + f72;
336
337	EvtComplex F78;
338	F78 = (-32.0/9.0)Lmu + 8.0EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0)
339	+ (-8.0EvtConst::pi/9.0)uniti +
340	(4.0/3.0EvtConst::piEvtConst::pi - 40.0/3.0)*shat +
341	(32.0EvtConst::piEvtConst::pi/9.0 - 316.0/9.0)shatshat +
342	(200.0EvtConst::piEvtConst::pi/27.0 - 658.0/9.0)shatshat*shat +
343	(-8.0logshat/9.0)(shat + shatshat + shatshat*shat);
344
345	c7eff = A7 - alphas/(4.0EvtConst::pi)(C1F71 + C2F72 + A8*F78);
346
347	return c7eff;
348	}
349
350
351	EvtComplex EvtbTosllAmp::GetC9Eff(double q2, bool nnlo, bool btod)
352	{
353
354	if (!nnlo) return 4.344;
355	double mbeff = 4.8;
356	double shat = q2/mbeff/mbeff;
357	double logshat;
358	logshat = log(shat);
359	double mchat = 0.29;
360
361
362	double muscale;
363	muscale = 2.5;
364	double alphas;
365	alphas = 0.267;
366	double A7;
367	A7 = -0.353 + 0.023;
368	double A8;
369	A8 = -0.164;
370	double A9;
371	A9 = 4.287 + (-0.218);
372	double A10;
373	A10 = -4.592 + 0.379;
374	double C1;
375	C1 = -0.697;
376	double C2;
377	C2 = 1.046;
378	double T9;
379	T9 = 0.114 + 0.280;
380	double U9;
381	U9 = 0.045 + 0.023;
382	double W9;
383	W9 = 0.044 + 0.016;
384
385	double Lmu;
386	Lmu = log(muscale/mbeff);
387
388
389	EvtComplex uniti(0.0,1.0);
390
391	EvtComplex hc;
392	double xarg;
393	xarg = 4.0*mchat/shat;
394	hc = -4.0/9.0log(mchatmchat) + 8.0/27.0 + 4.0*xarg/9.0;
395
396	if (xarg < 1.0)
397	{
398	hc = hc - 2.0/9.0(2.0 + xarg)sqrt(fabs(1.0 - xarg))*
399	(log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
400	uniti*EvtConst::pi);
401	}
402	else
403	{
404	hc = hc - 2.0/9.0(2.0 + xarg)sqrt(fabs(1.0 - xarg))*
405	2.0*atan(1.0/sqrt(xarg - 1.0));
406	}
407
408	EvtComplex h1;
409	xarg = 4.0/shat;
410	h1 = 8.0/27.0 + 4.0*xarg/9.0;
411	if (xarg < 1.0)
412	{
413	h1 = h1 - 2.0/9.0(2.0 + xarg)sqrt(fabs(1.0 - xarg))*
414	(log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
415	uniti*EvtConst::pi);
416	}
417	else
418	{
419	h1 = h1 - 2.0/9.0(2.0 + xarg)sqrt(fabs(1.0 - xarg))*
420	2.0*atan(1.0/sqrt(xarg - 1.0));
421	}
422
423
424	EvtComplex h0;
425	h0 = 8.0/27.0 - 4.0log(2.0)/9.0 + 4.0uniti*EvtConst::pi/9.0;
426
427
428	// X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
429	// h(\hat m_u^2, hat s))
430	EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
431	EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
432	EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
433	EvtComplex Vtb(1.0,0.0);
434
435	EvtComplex Xd;
436	Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0C1 + C2) (hc - h0);
437
438
439	EvtComplex c9eff=4.344;
440	if (shat > 0.25)
441	{
442	c9eff = A9 + T9hc + U9h1 + W9*h0;
443	if (btod)
444	{
445	c9eff += Xd;
446	}
447
448	return c9eff;
449	}
450
451	// change energy scale to 5.0 for full NNLO calculation below shat = 0.25
452	muscale = 5.0;
453	alphas = 0.215;
454	A9 = 4.174 + (-0.035);
455	C1 = -0.487;
456	C2 = 1.024;
457	A8 = -0.148;
458	T9 = 0.374 + 0.252;
459	U9 = 0.033 + 0.015;
460	W9 = 0.032 + 0.012;
461	Lmu = log(muscale/mbeff);
462
463	EvtComplex F91;
464	EvtComplex f91;
465	EvtComplex k9100(-11.973,0.16371);
466	EvtComplex k9101(-0.081271,-0.059691);
467	EvtComplex k9110(-28.432,-0.25044);
468	EvtComplex k9111(-0.040243,0.016442);
469	EvtComplex k9120(-57.114,-0.86486);
470	EvtComplex k9121(-0.035191,0.027909);
471	EvtComplex k9130(-128.8,-2.5243);
472	EvtComplex k9131(-0.017587,0.050639);
473	f91 = k9100 + k9101logshat + shat(k9110 + k9111*logshat) +
474	shatshat(k9120 + k9121*logshat) +
475	shatshatshat(k9130 + k9131logshat);
476	F91 = (-1424.0/729.0 + 16.0unitiEvtConst::pi/243.0
477	+ 64.0/27.0log(mchat))Lmu - 16.0Lmulogshat/243.0 +
478	(16.0/1215.0 - 32.0/135.0/mchat/mchat)Lmushat +
479	(4.0/2835.0 - 8.0/315.0/mchat/mchat/mchat/mchat)Lmushat*shat +
480	(16.0/76545.0 - 32.0/8505.0/mchat/mchat/mchat/mchat/mchat/mchat)*
481	Lmushatshatshat -256.0Lmu*Lmu/243.0 + f91;
482
483	EvtComplex F92;
484	EvtComplex f92;
485	EvtComplex k9200(6.6338,-0.98225);
486	EvtComplex k9201(0.48763,0.35815);
487	EvtComplex k9210(3.3585,1.5026);
488	EvtComplex k9211(0.24146,-0.098649);
489	EvtComplex k9220(-1.1906,5.1892);
490	EvtComplex k9221(0.21115,-0.16745);
491	EvtComplex k9230(-17.12,15.146);
492	EvtComplex k9231(0.10552,-0.30383);
493	f92 = k9200 + k9201logshat + shat(k9210 + k9211*logshat) +
494	shatshat(k9220 + k9221*logshat) +
495	shatshatshat(k9230 + k9231logshat);
496	F92 = (256.0/243.0 - 32.0unitiEvtConst::pi/81.0
497	- 128.0/9.0log(mchat))Lmu + 32.0Lmulogshat/81.0 +
498	(-32.0/405.0 + 64.0/45.0/mchat/mchat)Lmushat +
499	(-8.0/945.0 + 16.0/105.0/mchat/mchat/mchat/mchat)Lmushat*shat +
500	(-32.0/25515.0 + 64.0/2835.0/mchat/mchat/mchat/mchat/mchat/mchat)*
501	Lmushatshatshat + 512.0Lmu*Lmu/81.0 + f92;
502
503	EvtComplex F98;
504	F98 = 104.0/9.0 - 32.0EvtConst::piEvtConst::pi/27.0 +
505	(1184.0/27.0 - 40.0EvtConst::piEvtConst::pi/9.0)*shat +
506	(14212.0/135.0 - 32.0EvtConst::piEvtConst::pi/3.0)shatshat +
507	(193444.0/945.0 - 560.0EvtConst::piEvtConst::pi/27.0)shatshat*shat +
508	16.0logshat/9.0(1.0 + shat + shatshat + shatshat*shat);
509
510	Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0C1 + C2) (hc - h0);
511
512	c9eff = A9 + T9hc + U9h1 + W9*h0 -
513	alphas/(4.0EvtConst::pi)(C1F91 + C2F92 + A8*F98);
514	if (btod)
515	{
516	c9eff += Xd;
517	}
518
519	return c9eff;
520	}
521
522	EvtComplex EvtbTosllAmp::GetC10Eff(double /q2/, bool nnlo)
523	{
524
525	if (!nnlo) return -4.669;
526	double A10;
527	A10 = -4.592 + 0.379;
528
529	EvtComplex c10eff;
530	c10eff = A10;
531
532	return c10eff;
533	}
534
535	double EvtbTosllAmp::dGdsProb(double mb, double ms, double ml,
536	double s)
537	{
538	// Compute the decay probability density function given a value of s
539	// according to Ali's paper
540
541
542	double delta, lambda, prob;
543	double f1, f2, f3, f4;
544	double msh, mlh, sh;
545
546	mlh = ml / mb;
547	msh = ms / mb;
548	sh = s / (mb*mb);
549
550	EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
551	EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
552	EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);
553
554	double alphas = 0.119/
555	(1 + 0.119log(pow(4.8,2)/pow(91.1867,2))23.0/12.0/EvtConst::pi);
556	double omega9 = -2.0/9.0EvtConst::piEvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
557	- 2.0/3.0log(sh)log(1.0-sh)
558	- (5.0+4.0sh)/(3.0(1.0+2.0sh)) log(1.0-sh)
559	- 2.0sh(1.0+sh)(1.0-2.0sh)
560	/(3.0pow(1.0-sh,2)(1.0+2.0sh)) log(sh)
561	+ (5.0+9.0sh-6.0shsh)/(6.0(1.0-sh)(1.0+2.0sh));
562	double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
563	double omega7 = -8.0/3.0*log(4.8/mb)
564	-4.0/3.0*EvtDiLog::DiLog(sh)
565	-2.0/9.0EvtConst::piEvtConst::pi
566	-2.0/3.0log(sh)log(1.0-sh)
567	-log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
568	-2.0/3.0sh(2.0 - 2.0sh - shsh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
569	-(16.0 - 11.0sh - 17.0sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
570	double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
571
572	double omega79 = -4.0/3.0*log(4.8/mb)
573	-4.0/3.0*EvtDiLog::DiLog(sh)
574	-2.0/9.0EvtConst::piEvtConst::pi
575	-2.0/3.0log(sh)log(1.0-sh)
576	-1.0/9.0(2.0+7.0sh)*log(1.0 - sh)/sh
577	-2.0/9.0sh(3.0 - 2.0sh)log(sh)/pow((1.0 - sh),2)
578	+1.0/18.0(5.0 - 9.0sh)/(1.0 - sh);
579	double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
580
581	double c7c9 = abs(c7eff)*real(c9eff);
582	c7c9 *= pow(eta79,2);
583	double c7c7 = pow(abs(c7eff),2);
584	c7c7 *= pow(eta7,2);
585
586	double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
587	c9c9plusc10c10 *= pow(eta9,2);
588	double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
589	c9c9minusc10c10 *= pow(eta9,2);
590
591	lambda = 1.0 + shsh + pow(msh,4) - 2.0(sh + shmshmsh + msh*msh);
592
593	f1 = pow(1.0-mshmsh,2) - sh(1.0 + msh*msh);
594	f2 = 2.0(1.0 + mshmsh) * pow(1.0-msh*msh,2)
595	- sh(1.0 + 14.0mshmsh + pow(msh,4)) - shsh(1.0 + mshmsh);
596	f3 = pow(1.0-mshmsh,2) + sh(1.0 + mshmsh) - 2.0sh*sh
597	+ lambda2.0mlh*mlh/sh;
598	f4 = 1.0 - sh + msh*msh;
599
600	delta = ( 12.0c7c9f1 + 4.0c7c7f2/sh ) * (1.0 + 2.0mlhmlh/sh)
601	+ c9c9plusc10c10*f3
602	+ 6.0mlhmlhc9c9minusc10c10f4;
603
604	prob = sqrt(lambda(1.0 - 4.0mlhmlh/sh)) delta;
605
606	return prob;
607	}
608
609	double EvtbTosllAmp::dGdsdupProb(double mb, double ms, double ml,
610	double s, double u)
611	{
612	// Compute the decay probability density function given a value of s and u
613	// according to Ali's paper
614
615	double prob;
616	double f1sp, f2sp, f3sp;
617
618	double sh = s / (mb*mb);
619
620	EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
621	EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
622	EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);
623
624	double alphas = 0.119/
625	(1 + 0.119log(pow(4.8,2)/pow(91.1867,2))23.0/12.0/EvtConst::pi);
626	double omega9 = - 2.0/9.0EvtConst::piEvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
627	- 2.0/3.0log(sh)log(1.0-sh)
628	- (5.0+4.0sh)/(3.0(1.0+2.0sh)) log(1.0-sh)
629	- 2.0sh(1.0+sh)(1.0-2.0sh)
630	/(3.0pow(1.0-sh,2)(1.0+2.0sh)) log(sh)
631	+ (5.0+9.0sh-6.0shsh)/(6.0(1.0-sh)(1.0+2.0sh));
632	double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
633	double omega7 = -8.0/3.0*log(4.8/mb)
634	-4.0/3.0*EvtDiLog::DiLog(sh)
635	-2.0/9.0EvtConst::piEvtConst::pi
636	-2.0/3.0log(sh)log(1.0-sh)
637	-log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
638	-2.0/3.0sh(2.0 - 2.0sh - shsh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
639	-(16.0 - 11.0sh - 17.0sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
640	double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
641
642	double omega79 = -4.0/3.0*log(4.8/mb)
643	-4.0/3.0*EvtDiLog::DiLog(sh)
644	-2.0/9.0EvtConst::piEvtConst::pi
645	-2.0/3.0log(sh)log(1.0-sh)
646	-1.0/9.0(2.0+7.0sh)*log(1.0 - sh)/sh
647	-2.0/9.0sh(3.0 - 2.0sh)log(sh)/pow((1.0 - sh),2)
648	+1.0/18.0(5.0 - 9.0sh)/(1.0 - sh);
649	double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
650
651	double c7c9 = abs(c7eff)*real(c9eff);
652	c7c9 *= pow(eta79,2);
653	double c7c7 = pow(abs(c7eff),2);
654	c7c7 *= pow(eta7,2);
655
656	double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
657	c9c9plusc10c10 *= pow(eta9,2);
658	double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
659	c9c9minusc10c10 *= pow(eta9,2);
660	double c7c10 = abs(c7eff)*real(c10eff);
661	c7c10 = eta7; c7c10 = eta9;
662	double c9c10 = real(c9eff)*real(c10eff);
663	c9c10 *= pow(eta9,2);
664
665	f1sp = ( pow(mbmb-msms,2) - ss) c9c9plusc10c10
666	+ 4.0( pow(mb,4) - msmsmbmb - pow(ms,4)(1.0 - msms/(mb*mb))
667	- 8.0smsms - ss(1.0 + msms/(mbmb) ))mbmbc7c7/s
668	// kludged mass term
669	(1.0 + 2.0ml*ml/s)
670	- 8.0(s(mbmb + msms) - pow(mbmb-msms,2)) * c7c9
671	// kludged mass term
672	(1.0 + 2.0ml*ml/s);
673
674	f2sp = 4.0sc9c10 + 8.0(mbmb + msms)c7c10;
675	f3sp = - (c9c9plusc10c10)
676	+ 4.0(1.0 + pow(ms/mb,4)) mbmbc7c7/s
677	// kludged mass term
678	(1.0 + 2.0ml*ml/s);
679
680	prob = (f1sp + f2spu + f3spu*u)/ pow(mb,3);
681
682	return prob;
683	}
684
685
686
687
688
689
690
691
692
693
694
695
696