EvtSemiLeptonicAmp.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// 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: EvtSemiLeptonicAmp.cc
//
// Description: Routine to implement semileptonic decays to pseudo-scalar
//              mesons.
//
// Modification history:
//
//    DJL       April 17,1998       Module created
//
//------------------------------------------------------------------------
//
#include "EvtAmp.hh"
#include "EvtDiracSpinor.hh"
#include "EvtGenKine.hh"
#include "EvtId.hh"
#include "EvtPDL.hh"
#include "EvtParticle.hh"
#include "EvtPatches.hh"
#include "EvtReport.hh"
#include "EvtScalarParticle.hh"
#include "EvtSemiLeptonicScalarAmp.hh"
#include "EvtTensor4C.hh"
#include "EvtTensorParticle.hh"
#include "EvtVector4C.hh"
#include "EvtVectorParticle.hh"
 
double EvtSemiLeptonicAmp::CalcMaxProb( EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug,
                                        EvtSemiLeptonicFF* FormFactors ) {
 
  // 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->set_type(EvtSpinType::SCALAR);
  root_part->setDiagonalSpinDensity();
 
  EvtParticle *daughter, *lep, *trino;
 
  EvtAmp amp;
 
  EvtId listdaug[3];
  listdaug[0] = meson;
  listdaug[1] = lepton;
  listdaug[2] = nudaug;
 
  amp.init( parent, 3, listdaug );
 
  root_part->makeDaughters( 3, listdaug );
  daughter = root_part->getDaug( 0 );
  lep      = root_part->getDaug( 1 );
  trino    = root_part->getDaug( 2 );
 
  // cludge to avoid generating random numbers!
  daughter->noLifeTime();
  lep->noLifeTime();
  trino->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, p4lepton, p4nu, p4w;
  double      q2max;
 
  double q2, elepton, plepton;
  int    i, j;
  double erho, prho, costl;
 
  double maxfoundprob = 0.0;
  double prob         = -10.0;
  int    massiter;
 
  for ( massiter = 0; massiter < 3; massiter++ )
  {
 
    mass[0] = EvtPDL::getMeanMass( meson );
    mass[1] = EvtPDL::getMeanMass( lepton );
    mass[2] = EvtPDL::getMeanMass( nudaug );
    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
 
    for ( i = 0; i < 25; i++ )
    {
      q2 = ( ( i + 0.5 ) * q2max ) / 25.0;
 
      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.
        p4lepton.set( elepton, 0.0, plepton * sqrt( 1.0 - costl * costl ), plepton * costl );
        p4nu.set( plepton, 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 );
        p4lepton = boostTo( p4lepton, boost );
        p4nu     = boostTo( p4nu, boost );
 
        // Now initialize the daughters...
 
        daughter->init( meson, p4meson );
        lep->init( lepton, p4lepton );
        trino->init( nudaug, p4nu );
 
        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!!
        prob = rho.NormalizedProb( amp.getSpinDensity() );
 
        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 ( prob > maxfoundprob ) { maxfoundprob = prob; }
    }
    if ( EvtPDL::getWidth( meson ) <= 0.0 )
    {
      // if the particle is narrow dont bother with changing the mass.
      massiter = 4;
    }
  }
  root_part->deleteTree();
 
  maxfoundprob *= 1.1;
  return maxfoundprob;
}