EvtBTo4piCP.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: EvtBTo4piCP.cc
//
// Description: Routine to decay B->pi+ pi- pi+ pi-.
//
// Modification history:
//
//    RYD     March 2, 1997         Module created
//
//------------------------------------------------------------------------
//
#include "EvtBTo4piCP.hh"
#include "../EvtGenBase/EvtCPUtil.hh"
#include "../EvtGenBase/EvtConst.hh"
#include "../EvtGenBase/EvtGenKine.hh"
#include "../EvtGenBase/EvtId.hh"
#include "../EvtGenBase/EvtPDL.hh"
#include "../EvtGenBase/EvtParticle.hh"
#include "../EvtGenBase/EvtPatches.hh"
#include "../EvtGenBase/EvtReport.hh"
#include <stdlib.h>
#include <string>
 
EvtBTo4piCP::~EvtBTo4piCP() {}
 
EvtComplex EvtAmpA2( const EvtVector4R& p4pi1, const EvtVector4R& p4pi2,
                     const EvtVector4R& p4pi3, const EvtVector4R& p4pi4 ) {
 
  // added by Lange Jan4,2000
  static EvtId A2M  = EvtPDL::getId( "a_2-" );
  static EvtId RHO0 = EvtPDL::getId( "rho0" );
 
  EvtVector4R p4a2, p4rho, p4b;
 
  p4rho = p4pi1 + p4pi2;
 
  p4a2 = p4rho + p4pi3;
 
  p4b = p4a2 + p4pi4;
 
  EvtVector4R p4b_a2, p4rho_a2, p4pi1_a2, p4a2_a2;
 
  p4b_a2   = boostTo( p4b, p4a2 );
  p4rho_a2 = boostTo( p4rho, p4a2 );
  p4pi1_a2 = boostTo( p4pi1, p4a2 );
  p4a2_a2  = boostTo( p4a2, p4a2 );
 
  EvtVector4R p4pi1_rho;
 
  p4pi1_rho = boostTo( p4pi1_a2, p4rho_a2 );
 
  EvtVector4R vb, vrho, vpi, t;
 
  vb   = p4b_a2 / p4b_a2.d3mag();
  vrho = p4rho_a2 / p4rho_a2.d3mag();
  vpi  = p4pi1_rho / p4pi1_rho.d3mag();
 
  t.set( 1.0, 0.0, 0.0, 0.0 );
 
  //  EvtComplex amp_a1,amp_a2;
  EvtComplex amp_a2;
 
  //  double bwm_a1=EvtPDL::getMeanMass(A1M);
  //  double gamma_a1=EvtPDL::getWidth(A1M);
  double bwm_a2    = EvtPDL::getMeanMass( A2M );
  double gamma_a2  = EvtPDL::getWidth( A2M );
  double bwm_rho   = EvtPDL::getMeanMass( RHO0 );
  double gamma_rho = EvtPDL::getWidth( RHO0 );
 
  amp_a2 = ( sqrt( gamma_a2 / EvtConst::twoPi ) /
             ( ( p4a2 ).mass() - bwm_a2 - EvtComplex( 0.0, 0.5 * gamma_a2 ) ) ) *
           ( sqrt( gamma_rho / EvtConst::twoPi ) /
             ( ( p4rho ).mass() - bwm_rho - EvtComplex( 0.0, 0.5 * gamma_rho ) ) );
 
  return amp_a2 *
         ( vb.get( 1 ) * vrho.get( 1 ) + vb.get( 2 ) * vrho.get( 2 ) +
           vb.get( 3 ) * vrho.get( 3 ) ) *
         ( vpi.get( 1 ) * ( vb.get( 2 ) * vrho.get( 3 ) - vb.get( 3 ) * vrho.get( 2 ) ) +
           vpi.get( 2 ) * ( vb.get( 3 ) * vrho.get( 1 ) - vb.get( 1 ) * vrho.get( 3 ) ) +
           vpi.get( 3 ) * ( vb.get( 1 ) * vrho.get( 2 ) - vb.get( 2 ) * vrho.get( 1 ) ) );
}
 
EvtComplex EvtAmpA1( const EvtVector4R& p4pi1, const EvtVector4R& p4pi2,
                     const EvtVector4R& p4pi3, const EvtVector4R& p4pi4 ) {
 
  // added by Lange Jan4,2000
  static EvtId A1M  = EvtPDL::getId( "a_1-" );
  static EvtId RHO0 = EvtPDL::getId( "rho0" );
 
  EvtVector4R p4a1, p4rho, p4b;
 
  p4rho = p4pi1 + p4pi2;
 
  p4a1 = p4rho + p4pi3;
 
  p4b = p4a1 + p4pi4;
 
  EvtVector4R p4b_a1, p4rho_a1, p4pi1_a1, p4a1_a1;
 
  p4b_a1   = boostTo( p4b, p4a1 );
  p4rho_a1 = boostTo( p4rho, p4a1 );
  p4pi1_a1 = boostTo( p4pi1, p4a1 );
  p4a1_a1  = boostTo( p4a1, p4a1 );
 
  EvtVector4R p4pi1_rho;
 
  p4pi1_rho = boostTo( p4pi1_a1, p4rho_a1 );
 
  EvtVector4R vb, vrho, vpi, t;
 
  vb   = p4b_a1 / p4b_a1.d3mag();
  vrho = p4rho_a1 / p4rho_a1.d3mag();
  vpi  = p4pi1_rho / p4pi1_rho.d3mag();
 
  t.set( 1.0, 0.0, 0.0, 0.0 );
 
  EvtComplex amp_a1;
 
  double bwm_a1   = EvtPDL::getMeanMass( A1M );
  double gamma_a1 = EvtPDL::getWidth( A1M );
  //  double bwm_a2=EvtPDL::getMeanMass(A2M);
  //  double gamma_a2=EvtPDL::getWidth(A2M);
  double bwm_rho   = EvtPDL::getMeanMass( RHO0 );
  double gamma_rho = EvtPDL::getWidth( RHO0 );
 
  amp_a1 = ( sqrt( gamma_a1 / EvtConst::twoPi ) /
             ( ( p4a1 ).mass() - bwm_a1 - EvtComplex( 0.0, 0.5 * gamma_a1 ) ) ) *
           ( sqrt( gamma_rho / EvtConst::twoPi ) /
             ( ( p4rho ).mass() - bwm_rho - EvtComplex( 0.0, 0.5 * gamma_rho ) ) );
 
  return amp_a1 * ( vb.get( 1 ) * vpi.get( 1 ) + vb.get( 2 ) * vpi.get( 2 ) +
                    vb.get( 3 ) * vpi.get( 3 ) );
}
 
void EvtBTo4piCP::getName( std::string& model_name ) { model_name = "BTO4PI_CP"; }
 
EvtDecayBase* EvtBTo4piCP::clone() { return new EvtBTo4piCP; }
 
void EvtBTo4piCP::init() {
 
  // check that there are 18 arguments
  checkNArg( 18 );
  checkNDaug( 4 );
 
  checkSpinParent( EvtSpinType::SCALAR );
 
  checkSpinDaughter( 0, EvtSpinType::SCALAR );
  checkSpinDaughter( 1, EvtSpinType::SCALAR );
  checkSpinDaughter( 2, EvtSpinType::SCALAR );
  checkSpinDaughter( 3, EvtSpinType::SCALAR );
}
 
void EvtBTo4piCP::decay( EvtParticle* p ) {
 
  // added by Lange Jan4,2000
  static EvtId B0  = EvtPDL::getId( "B0" );
  static EvtId B0B = EvtPDL::getId( "anti-B0" );
 
  double t;
  EvtId  other_b;
 
  EvtCPUtil::OtherB( p, t, other_b );
 
  p->initializePhaseSpace( getNDaug(), getDaugs() );
  EvtVector4R mom1 = p->getDaug( 0 )->getP4();
  EvtVector4R mom2 = p->getDaug( 1 )->getP4();
  EvtVector4R mom3 = p->getDaug( 2 )->getP4();
  EvtVector4R mom4 = p->getDaug( 3 )->getP4();
 
  //  double alpha=getArg(0);
  // double dm=getArg(1);
 
  EvtComplex amp;
 
  EvtComplex A, Abar;
 
  EvtComplex A_a1p, Abar_a1p, A_a2p, Abar_a2p;
  EvtComplex A_a1m, Abar_a1m, A_a2m, Abar_a2m;
 
  A_a1p    = EvtComplex( getArg( 2 ) * cos( getArg( 3 ) ), getArg( 2 ) * sin( getArg( 3 ) ) );
  Abar_a1p = EvtComplex( getArg( 4 ) * cos( getArg( 5 ) ), getArg( 4 ) * sin( getArg( 5 ) ) );
 
  A_a2p    = EvtComplex( getArg( 6 ) * cos( getArg( 7 ) ), getArg( 6 ) * sin( getArg( 7 ) ) );
  Abar_a2p = EvtComplex( getArg( 8 ) * cos( getArg( 9 ) ), getArg( 8 ) * sin( getArg( 9 ) ) );
 
  A_a1m = EvtComplex( getArg( 10 ) * cos( getArg( 11 ) ), getArg( 10 ) * sin( getArg( 11 ) ) );
  Abar_a1m =
      EvtComplex( getArg( 12 ) * cos( getArg( 13 ) ), getArg( 12 ) * sin( getArg( 13 ) ) );
 
  A_a2m = EvtComplex( getArg( 14 ) * cos( getArg( 15 ) ), getArg( 14 ) * sin( getArg( 15 ) ) );
  Abar_a2m =
      EvtComplex( getArg( 16 ) * cos( getArg( 17 ) ), getArg( 16 ) * sin( getArg( 17 ) ) );
 
  EvtComplex a2p_amp = EvtAmpA2( mom1, mom2, mom3, mom4 ) +
                       EvtAmpA2( mom1, mom4, mom3, mom2 ) +
                       EvtAmpA2( mom3, mom2, mom1, mom4 ) + EvtAmpA2( mom3, mom4, mom1, mom2 );
 
  EvtComplex a2m_amp = EvtAmpA2( mom2, mom3, mom4, mom1 ) +
                       EvtAmpA2( mom2, mom1, mom4, mom3 ) +
                       EvtAmpA2( mom4, mom3, mom2, mom1 ) + EvtAmpA2( mom4, mom1, mom2, mom3 );
 
  EvtComplex a1p_amp = EvtAmpA1( mom1, mom2, mom3, mom4 ) +
                       EvtAmpA1( mom1, mom4, mom3, mom2 ) +
                       EvtAmpA1( mom3, mom2, mom1, mom4 ) + EvtAmpA1( mom3, mom4, mom1, mom2 );
 
  EvtComplex a1m_amp = EvtAmpA1( mom2, mom3, mom4, mom1 ) +
                       EvtAmpA1( mom2, mom1, mom4, mom3 ) +
                       EvtAmpA1( mom4, mom3, mom2, mom1 ) + EvtAmpA1( mom4, mom1, mom2, mom3 );
 
  A    = A_a2p * a2p_amp + A_a1p * a1p_amp + A_a2m * a2m_amp + A_a1m * a1m_amp;
  Abar = Abar_a2p * a2p_amp + Abar_a1p * a1p_amp + Abar_a2m * a2m_amp + Abar_a1m * a1m_amp;
 
  if ( other_b == B0B )
  {
    amp = A * cos( getArg( 1 ) * t / ( 2 * EvtConst::c ) ) +
          EvtComplex( cos( -2.0 * getArg( 0 ) ), sin( -2.0 * getArg( 0 ) ) ) * getArg( 2 ) *
              EvtComplex( 0.0, 1.0 ) * Abar * sin( getArg( 1 ) * t / ( 2 * EvtConst::c ) );
  }
  if ( other_b == B0 )
  {
    amp = A * EvtComplex( cos( 2.0 * getArg( 0 ) ), sin( 2.0 * getArg( 0 ) ) ) *
              EvtComplex( 0.0, 1.0 ) * sin( getArg( 1 ) * t / ( 2 * EvtConst::c ) ) +
          getArg( 2 ) * Abar * cos( getArg( 1 ) * t / ( 2 * EvtConst::c ) );
  }
 
  vertex( amp );
 
  return;
}