EvtSVSCPiso.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: EvtSVSCPiso.cc
//
// Description: Routine to decay scalar -> vectors scalar
//              with CP violation and isospin amplitudes.
//              More specifically, it is indended to handle
//              decays like B->rho pi and B->a1 pi.
//
// Modification history:
//
//    RYD/NK    Febuary 16, 1998          Module created
//
//------------------------------------------------------------------------
//
#include "EvtSVSCPiso.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/EvtRandom.hh"
#include "../EvtGenBase/EvtReport.hh"
#include "../EvtGenBase/EvtVector4C.hh"
#include <stdlib.h>
#include <string>
 
EvtSVSCPiso::~EvtSVSCPiso() {}
 
void EvtSVSCPiso::getName( std::string& model_name ) { model_name = "SVS_CP_ISO"; }
 
EvtDecayBase* EvtSVSCPiso::clone() { return new EvtSVSCPiso; }
 
void EvtSVSCPiso::init() {
 
  // check that there are 26 arguments
  checkNArg( 26 );
  checkNDaug( 2 );
 
  checkSpinParent( EvtSpinType::SCALAR );
 
  checkSpinDaughter( 0, EvtSpinType::VECTOR );
  checkSpinDaughter( 1, EvtSpinType::SCALAR );
}
 
void EvtSVSCPiso::initProbMax() {
 
  // this might need some revision..
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) > 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) == 0 ) )
  { setProbMax( 2.0 * ( getArg( 3 ) * getArg( 3 ) + 4.0 * getArg( 23 ) * getArg( 23 ) ) ); }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) < 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) == 0 ) )
  { setProbMax( 2.0 * ( getArg( 5 ) * getArg( 5 ) + 4.0 * getArg( 25 ) * getArg( 25 ) ) ); }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) == 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) > 0 ) )
  { setProbMax( 2.0 * ( getArg( 7 ) * getArg( 7 ) + 4.0 * getArg( 23 ) * getArg( 23 ) ) ); }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) == 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) < 0 ) )
  { setProbMax( 2.0 * ( getArg( 9 ) * getArg( 9 ) + 4.0 * getArg( 25 ) * getArg( 25 ) ) ); }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) > 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) < 0 ) )
  {
    setProbMax( 2.0 * ( getArg( 11 ) * getArg( 11 ) + getArg( 23 ) * getArg( 23 ) +
                        getArg( 19 ) * getArg( 19 ) + getArg( 13 ) * getArg( 13 ) +
                        getArg( 25 ) * getArg( 25 ) + getArg( 21 ) * getArg( 21 ) ) );
  }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) < 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) > 0 ) )
  {
    setProbMax( 2.0 * ( getArg( 15 ) * getArg( 15 ) + getArg( 23 ) * getArg( 23 ) +
                        getArg( 19 ) * getArg( 19 ) + getArg( 17 ) * getArg( 17 ) +
                        getArg( 25 ) * getArg( 25 ) + getArg( 21 ) * getArg( 21 ) ) );
  }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) == 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) == 0 ) )
  {
    setProbMax( 2.0 * ( getArg( 7 ) * getArg( 7 ) + getArg( 3 ) * getArg( 3 ) +
                        getArg( 11 ) * getArg( 11 ) + getArg( 15 ) * getArg( 15 ) +
                        4.0 * getArg( 19 ) * getArg( 19 ) + getArg( 9 ) * getArg( 9 ) +
                        getArg( 5 ) * getArg( 5 ) + getArg( 13 ) * getArg( 13 ) +
                        getArg( 17 ) * getArg( 17 ) + 4.0 * getArg( 21 ) * getArg( 21 ) ) );
  }
}
 
void EvtSVSCPiso::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;
  int    charged( 0 );
 
  int   first_time = 0;
  int   flip       = 0;
  EvtId ds[2];
 
  // randomly generate the tag (B0 or B0B)
 
  double tag = EvtRandom::Flat( 0.0, 1.0 );
  if ( tag < 0.5 )
  {
 
    EvtCPUtil::OtherB( p, t, other_b, 1.0 );
    other_b = B0;
  }
  else
  {
 
    EvtCPUtil::OtherB( p, t, other_b, 0.0 );
    other_b = B0B;
  }
 
  if ( p->getNDaug() == 0 ) first_time = 1;
 
  if ( first_time )
  {
    if ( EvtRandom::Flat( 0.0, 1.0 ) < getArg( 3 ) ) flip = 1;
  }
  else
  {
    if ( getDaug( 0 ) != p->getDaug( 0 )->getId() ) flip = 1;
  }
 
  if ( !flip )
  {
    ds[0] = getDaug( 0 );
    ds[1] = getDaug( 1 );
  }
  else
  {
    ds[0] = EvtPDL::chargeConj( getDaug( 0 ) );
    ds[1] = EvtPDL::chargeConj( getDaug( 1 ) );
  }
 
  p->initializePhaseSpace( getNDaug(), ds );
 
  EvtParticle *v, *s;
  v = p->getDaug( 0 );
  s = p->getDaug( 1 );
 
  EvtComplex amp;
 
  EvtComplex A_f, Abar_f;
  EvtComplex A_fbar, Abar_fbar;
  EvtComplex Apm, Apm_bar, Amp, Amp_bar;
 
  EvtComplex Tp0, Tp0_bar, T0p, T0p_bar, Tpm, Tpm_bar, Tmp, Tmp_bar;
  EvtComplex P1, P1_bar, P0, P0_bar;
 
  Tp0     = EvtComplex( getArg( 3 ) * cos( getArg( 4 ) ), getArg( 3 ) * sin( getArg( 4 ) ) );
  Tp0_bar = EvtComplex( getArg( 5 ) * cos( getArg( 6 ) ), getArg( 5 ) * sin( getArg( 6 ) ) );
  T0p     = EvtComplex( getArg( 7 ) * cos( getArg( 8 ) ), getArg( 7 ) * sin( getArg( 8 ) ) );
  T0p_bar = EvtComplex( getArg( 9 ) * cos( getArg( 10 ) ), getArg( 9 ) * sin( getArg( 10 ) ) );
  Tpm = EvtComplex( getArg( 11 ) * cos( getArg( 12 ) ), getArg( 11 ) * sin( getArg( 12 ) ) );
  Tpm_bar =
      EvtComplex( getArg( 13 ) * cos( getArg( 14 ) ), getArg( 13 ) * sin( getArg( 14 ) ) );
  Tmp = EvtComplex( getArg( 15 ) * cos( getArg( 16 ) ), getArg( 15 ) * sin( getArg( 16 ) ) );
  Tmp_bar =
      EvtComplex( getArg( 17 ) * cos( getArg( 18 ) ), getArg( 17 ) * sin( getArg( 18 ) ) );
  P0 = EvtComplex( getArg( 19 ) * cos( getArg( 20 ) ), getArg( 19 ) * sin( getArg( 20 ) ) );
  P0_bar =
      EvtComplex( getArg( 21 ) * cos( getArg( 22 ) ), getArg( 21 ) * sin( getArg( 22 ) ) );
  P1 = EvtComplex( getArg( 23 ) * cos( getArg( 24 ) ), getArg( 23 ) * sin( getArg( 24 ) ) );
  P1_bar =
      EvtComplex( getArg( 25 ) * cos( getArg( 26 ) ), getArg( 25 ) * sin( getArg( 26 ) ) );
 
  //***********************charged modes****************************
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) > 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) == 0 ) )
  {
 
    // V+ S0, so T+0 + 2 P1
 
    charged = 1;
    A_f     = Tp0 + 2.0 * P1;
  }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) < 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) == 0 ) )
  {
 
    // V- S0, so T+0_bar + 2P1_bar
 
    charged = 1;
    A_f     = Tp0_bar + 2.0 * P1_bar;
  }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) == 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) > 0 ) )
  {
 
    // V0 S+, so T0+ - 2 P1
 
    charged = 1;
    A_f     = T0p - 2.0 * P1;
  }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) == 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) < 0 ) )
  {
 
    // V0 S-, so T0+_bar - 2 P1_bar
 
    charged = 1;
    A_f     = T0p_bar - 2.0 * P1_bar;
  }
 
  //***********************neutral modes***************************
 
  // V+ S-, so Af = T+- + P1 + P0
  Apm     = Tpm + P1 + P0;
  Apm_bar = Tpm_bar + P1_bar + P0_bar;
 
  // V- S+, so Af = T-+ - P1 + P0
  Amp     = Tmp - P1 + P0;
  Amp_bar = Tmp_bar - P1_bar + P0;
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) > 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) < 0 ) )
  {
 
    // V+ S-
    charged   = 0;
    A_f       = Apm;
    Abar_f    = Apm_bar;
    A_fbar    = Amp;
    Abar_fbar = Amp_bar;
  }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) < 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) > 0 ) )
  {
 
    // V- S+
    charged   = 0;
    A_f       = Amp;
    Abar_f    = Amp_bar;
    A_fbar    = Apm;
    Abar_fbar = Apm_bar;
  }
 
  if ( ( EvtPDL::chg3( getDaug( 0 ) ) == 0 ) && ( EvtPDL::chg3( getDaug( 1 ) ) == 0 ) )
  {
 
    // V0 S0
    charged   = 0;
    A_f       = T0p + Tp0 - Tpm - Tmp - 2.0 * P0;
    Abar_f    = T0p_bar + Tp0_bar - Tpm_bar - Tmp_bar - 2.0 * P0_bar;
    A_fbar    = A_f;
    Abar_fbar = Abar_f;
  }
 
  if ( charged == 0 )
  {
 
    if ( !flip )
    {
      if ( other_b == B0B )
      {
 
        amp = A_f * cos( getArg( 1 ) * t / ( 2 * EvtConst::c ) ) +
              EvtComplex( cos( -2.0 * getArg( 0 ) ), sin( -2.0 * getArg( 0 ) ) ) *
                  EvtComplex( 0.0, 1.0 ) * Abar_f *
                  sin( getArg( 1 ) * t / ( 2 * EvtConst::c ) );
      }
      if ( other_b == B0 )
      {
 
        amp = A_f * EvtComplex( cos( 2.0 * getArg( 0 ) ), sin( 2.0 * getArg( 0 ) ) ) *
                  EvtComplex( 0.0, 1.0 ) * sin( getArg( 1 ) * t / ( 2 * EvtConst::c ) ) +
              Abar_f * cos( getArg( 1 ) * t / ( 2 * EvtConst::c ) );
      }
    }
    else
    {
      if ( other_b == B0B )
      {
 
        amp = A_fbar * cos( getArg( 1 ) * t / ( 2 * EvtConst::c ) ) +
              EvtComplex( cos( -2.0 * getArg( 0 ) ), sin( -2.0 * getArg( 0 ) ) ) *
                  EvtComplex( 0.0, 1.0 ) * Abar_fbar *
                  sin( getArg( 1 ) * t / ( 2 * EvtConst::c ) );
      }
      if ( other_b == B0 )
      {
 
        amp = A_fbar * EvtComplex( cos( 2.0 * getArg( 0 ) ), sin( 2.0 * getArg( 0 ) ) ) *
                  EvtComplex( 0.0, 1.0 ) * sin( getArg( 1 ) * t / ( 2 * EvtConst::c ) ) +
              Abar_fbar * cos( getArg( 1 ) * t / ( 2 * EvtConst::c ) );
      }
    }
  }
  else amp = A_f;
 
  EvtVector4R p4_parent;
 
  p4_parent = v->getP4() + s->getP4();
 
  double norm = 1.0 / v->getP4().d3mag();
 
  vertex( 0, amp * norm * p4_parent * ( v->epsParent( 0 ) ) );
  vertex( 1, amp * norm * p4_parent * ( v->epsParent( 1 ) ) );
  vertex( 2, amp * norm * p4_parent * ( v->epsParent( 2 ) ) );
 
  return;
}