EvtDsTophienu.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
// Environment:
//      This software is part of models developed at BES collaboration
//      based on the EvtGen framework.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/BesCopyright
//      Copyright (A) 2006      Ping Rong-Gang @IHEP
//
// Module: EvtDsTophienu.cc
//      the necessary file: EvtDsTophienu.hh
//
// Description: Ds+ -> phi e+ nu, phi -> pi+ pi- pi0
//
// Modification history:
//      Liaoyuan Dong    May 29 00:41:03 2023   Module created
//------------------------------------------------------------------------
#include "EvtDsTophienu.hh"
#include "../EvtGenBase/EvtComplex.hh"
#include "../EvtGenBase/EvtDecayTable.hh"
#include "../EvtGenBase/EvtGenKine.hh"
#include "../EvtGenBase/EvtPDL.hh"
#include "../EvtGenBase/EvtParticle.hh"
#include "../EvtGenBase/EvtPatches.hh"
#include "../EvtGenBase/EvtReport.hh"
#include "../EvtGenBase/EvtVector3R.hh"
#include <stdlib.h>
 
EvtDsTophienu::~EvtDsTophienu() {}
 
void EvtDsTophienu::getName( std::string& model_name ) { model_name = "DsTophienu"; }
 
EvtDecayBase* EvtDsTophienu::clone() { return new EvtDsTophienu; }
 
void EvtDsTophienu::init() {
  checkNArg( 0 );
  checkNDaug( 5 ); //  Ds+ -> pi+ pi- pi0 e+ nu
                   //  0   -> 1   2   3   4  5
  checkSpinParent( EvtSpinType::SCALAR );
  checkSpinDaughter( 0, EvtSpinType::SCALAR );
  checkSpinDaughter( 1, EvtSpinType::SCALAR );
  checkSpinDaughter( 2, EvtSpinType::SCALAR );
 
  std::cout << "Initializing EvtDsTophienu" << std::endl;
  mV   = 1.81;
  mA   = 2.61;
  V_0  = 1.411;
  A1_0 = 1;
  A2_0 = 0.788;
  mV2  = mV * mV;
  mA2  = mA * mA;
 
  m_phi  = 1.019461; // phi
  w_phi  = 0.004249;
  m_rho  = 0.77526; // rho
  w_rho  = 0.14910;
  rBW    = 3.0;
  rBW2   = rBW * rBW;
  mw_phi = m_phi * w_phi;
  m2_phi = m_phi * m_phi;
  m2_rho = m_rho * m_rho;
 
  mDs   = 1.9683;
  m2Ds  = mDs * mDs;
  mpi   = 0.13957;
  m2pi  = mpi * mpi;
  mpi0  = 0.1349766;
  m2pi0 = mpi0 * mpi0;
 
  Pi      = atan2( 0.0, -1.0 );
  root2   = sqrt( 2. );
  root2d3 = sqrt( 2. / 3 );
  root1d2 = sqrt( 0.5 );
  root3d2 = sqrt( 1.5 );
}
 
void EvtDsTophienu::initProbMax() { setProbMax( 240000.0 ); }
 
void EvtDsTophienu::decay( EvtParticle* p ) {
  /*
     double maxprob = 0.0;
     for(int ir=0;ir<=60000000;ir++){
        p->initializePhaseSpace(getNDaug(),getDaugs());
        EvtVector4R _pip = p->getDaug(0)->getP4(); // pi+
        EvtVector4R _pim = p->getDaug(1)->getP4(); // pi-
        EvtVector4R _pi0 = p->getDaug(2)->getP4(); // pi0
        EvtVector4R _e   = p->getDaug(3)->getP4(); // e+
        EvtVector4R _nu  = p->getDaug(4)->getP4(); // nu
        int pid = EvtPDL::getStdHep(p->getDaug(3)->getId());
        int charm;
        if(pid == -11) charm = 1;
        else charm = -1;
        double m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V;
 
        KinVGen(_pip, _pim, _pi0, _e, _nu, charm, m2, q2, cosV, cosL, chi, s12, s13, s23,
     spin_V); double _prob = calPDF(m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V);
        if(_prob>=maxprob) {
            maxprob=_prob;
            std::cout << ir << " prob= " << _prob << std::endl;
            cout << "EvtVector4R pip" << _pip << ";"<< endl;
            cout << "EvtVector4R pim" << _pim << ";"<< endl;
            cout << "EvtVector4R pi0" << _pi0 << ";"<< endl;
            cout << "EvtVector4R e"   << _e   << ";"<< endl;
            cout << "EvtVector4R nu"  << _nu  << ";"<< endl;
        }
     }
     std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
  */
  p->initializePhaseSpace( getNDaug(), getDaugs() );
  EvtVector4R pip = p->getDaug( 0 )->getP4();
  EvtVector4R pim = p->getDaug( 1 )->getP4();
  EvtVector4R pi0 = p->getDaug( 2 )->getP4();
  EvtVector4R e   = p->getDaug( 3 )->getP4();
  EvtVector4R nu  = p->getDaug( 4 )->getP4();
 
  int pid = EvtPDL::getStdHep( p->getDaug( 3 )->getId() );
  int charm;
  if ( pid == -11 ) charm = 1;
  else charm = -1;
  double m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V;
  KinVGen( pip, pim, pi0, e, nu, charm, m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V );
  double prob = calPDF( m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V );
  setProb( prob );
  return;
}
 
void EvtDsTophienu::KinVGen( EvtVector4R vp4_pip, EvtVector4R vp4_pim, EvtVector4R vp4_pi0,
                             EvtVector4R vp4_Lep, EvtVector4R vp4_Nu, int charm, double& m2,
                             double& q2, double& cosV, double& cosL, double& chi, double& s12,
                             double& s13, double& s23, double& spin_V ) {
  EvtVector4R vp4_3pi  = vp4_pip + vp4_pim + vp4_pi0;
  EvtVector4R vp4_rho0 = vp4_pip + vp4_pim;
  EvtVector4R vp4_rhop = vp4_pip + vp4_pi0;
  EvtVector4R vp4_rhom = vp4_pim + vp4_pi0;
  EvtVector4R vp4_W    = vp4_Lep + vp4_Nu;
 
  m2  = vp4_3pi.mass2();
  q2  = vp4_W.mass2();
  s12 = vp4_rho0.mass2();
  s13 = vp4_rhop.mass2();
  s23 = vp4_rhom.mass2();
 
  EvtVector4R boost;
  boost.set( vp4_3pi.get( 0 ), -vp4_3pi.get( 1 ), -vp4_3pi.get( 2 ), -vp4_3pi.get( 3 ) );
  EvtVector4R vp4_rhob = boostTo( vp4_rho0, boost );
  cosV                 = vp4_rhob.dot( vp4_3pi ) / ( vp4_rhob.d3mag() * vp4_3pi.d3mag() );
 
  EvtVector4R vp4_pipb = boostTo( vp4_pip, boost );
  EvtVector4R vp4_pimb = boostTo( vp4_pim, boost );
  EvtVector4R R        = vp4_pipb.cross( vp4_pimb );
  spin_V               = R.d3mag(); // sart(|pi+ times pi-|^2)
 
  boost.set( vp4_W.get( 0 ), -vp4_W.get( 1 ), -vp4_W.get( 2 ), -vp4_W.get( 3 ) );
  EvtVector4R vp4_Lepp = boostTo( vp4_Lep, boost );
  cosL                 = vp4_Lepp.dot( vp4_W ) / ( vp4_Lepp.d3mag() * vp4_W.d3mag() );
 
  EvtVector4R V = vp4_3pi / vp4_3pi.d3mag();
  EvtVector4R C = vp4_rhob.cross( V );
  C /= C.d3mag();
  EvtVector4R D = vp4_Lepp.cross( V );
  D /= D.d3mag();
  double sinx = C.cross( V ).dot( D );
  double cosx = C.dot( D );
  chi         = sinx > 0 ? acos( cosx ) : -acos( cosx );
  if ( charm == -1 ) chi = -chi;
}
 
double EvtDsTophienu::calPDF( double m2, double q2, double cosV, double cosL, double chi,
                              double s12, double s13, double s23, double spin_V ) {
  double m   = sqrt( m2 );
  double q   = sqrt( q2 );
  double m12 = sqrt( s12 ); // m(pi+pi-)
  double m13 = sqrt( s13 ); // m(pi+pi0)
  double m23 = sqrt( s23 ); // m(pi-pi0)
  // cout << "m12 = " << m12 << " m13 = " << m13 << " m23 = " << m23 << " m = " << m << " q = "
  // << q << endl;
 
  EvtComplex A_rhopi( 0.0, 0.0 );
  ResonanceV( s12, m12, s13, m13, s23, m23, A_rhopi );
  double A2_phi = spin_V * spin_V * abs2( A_rhopi ); // Eq.(1)
  // cout << "spin_V = " << spin_V << "2 = " << spin_V*spin_V << " A_rho = " << abs2(A_rhopi)
  // << endl;
 
  EvtComplex f11( 0.0, 0.0 );
  EvtComplex f21( 0.0, 0.0 );
  EvtComplex f31( 0.0, 0.0 );
  ResonanceP( m2, m, q2, q, s12, m12, s13, m13, s23, m23, f11, f21, f31 );
 
  double cosV2 = cosV * cosV;
  double sinV  = sqrt( 1.0 - cosV2 );
  double sinV2 = sinV * sinV;
 
  EvtComplex F1 = f11 * cosV;
  EvtComplex F2 = f21 * root1d2;
  EvtComplex F3 = f31 * root1d2;
  // cout << "F1= " << F1 << abs2(F1) << endl;
  // cout << "F2= " << F2 << abs2(F2) << endl;
  // cout << "F3= " << F3 << abs2(F3) << endl;
 
  double I1 = 0.25 * ( abs2( F1 ) + 1.5 * sinV2 * ( abs2( F2 ) + abs2( F3 ) ) );
  double I2 = -0.25 * ( abs2( F1 ) - 0.5 * sinV2 * ( abs2( F2 ) + abs2( F3 ) ) );
  double I3 = -0.25 * ( abs2( F2 ) - abs2( F3 ) ) * sinV2;
  double I4 = real( conj( F1 ) * F2 ) * sinV * 0.5;
  double I5 = real( conj( F1 ) * F3 ) * sinV;
  double I6 = real( conj( F2 ) * F3 ) * sinV2;
  // cout << "Ix = " << I1 << ", " << I2 << ", " << I3 << ", " << I4 << ", " << I5 << ", " <<
  // I6 << endl;
 
  double coschi  = cos( chi );
  double sinchi  = sin( chi );
  double sin2chi = 2.0 * sinchi * coschi;
  double cos2chi = 1.0 - 2.0 * sinchi * sinchi;
 
  double sinL  = sqrt( 1. - cosL * cosL );
  double sinL2 = sinL * sinL;
  double sin2L = 2.0 * sinL * cosL;
  double cos2L = 1.0 - 2.0 * sinL2;
 
  double I = A2_phi * ( I1 + I2 * cos2L + I3 * sinL2 * cos2chi + I4 * sin2L * coschi +
                        I5 * sinL * coschi + I6 * cosL );
  return I;
}
 
void EvtDsTophienu::ResonanceV( double s12, double m12, double s13, double m13, double s23,
                                double m23, EvtComplex& A_rhopi ) {
  double G_rho0 =
      w_rho / m12 *
      pow( ( s12 - 4.0 * m2pi0 ) / ( m2_rho - 4.0 * m2pi0 ), 1.5 ); // Gamma(rho0) Eq.(10)
  double G_rhop =
      w_rho / m13 *
      pow( ( s13 - 4.0 * m2pi ) / ( m2_rho - 4.0 * m2pi ), 1.5 ); // Gamma(rho+) Eq.(10)
  double G_rhom =
      w_rho / m23 *
      pow( ( s23 - 4.0 * m2pi ) / ( m2_rho - 4.0 * m2pi ), 1.5 ); // Gamma(rho-) Eq.(10)
 
  double     AA0  = s12 / m2_rho - 1.0;
  EvtComplex amp0 = EvtComplex( 1.0, 0.0 ) / EvtComplex( AA0, G_rho0 );
  double     AAp  = s13 / m2_rho - 1.0;
  EvtComplex ampp = EvtComplex( 1.0, 0.0 ) / EvtComplex( AAp, G_rhop );
  double     AAm  = s23 / m2_rho - 1.0;
  EvtComplex ampm = EvtComplex( 1.0, 0.0 ) / EvtComplex( AAm, G_rhom );
 
  A_rhopi = amp0 + ampp + ampm; // Eq.(10)
}
 
void EvtDsTophienu::ResonanceP( double m2, double m, double q2, double q, double s12,
                                double m12, double s13, double m13, double s23, double m23,
                                EvtComplex& F11, EvtComplex& F21, EvtComplex& F31 ) {
  double p_phi   = getPStar( m2Ds, m2, q2 );
  double summDsm = mDs + m;
  double V       = V_0 / ( 1.0 - q2 / ( mV2 ) );
  double A1      = A1_0 / ( 1.0 - q2 / ( mA2 ) );
  double A2      = A2_0 / ( 1.0 - q2 / ( mA2 ) );
  double A       = summDsm * A1;
  double B       = 2.0 * mDs * p_phi / summDsm * V;
 
  // construct the helicity form factor
  double H0 =
      0.5 / ( m * q ) *
      ( ( m2Ds - m2 - q2 ) * summDsm * A1 - 4.0 * ( m2Ds * p_phi * p_phi ) / summDsm * A2 );
  double Hp = A - B;
  double Hm = A + B;
 
  // calculate alpha
  double pStar0 = getPStar( m2_phi, s12, m2pi0 );
  // double pStar0 = getPStar(m2_phi, m2_rho, m2pi0);
  double alpha = sqrt( 3. * Pi * 0.154 / ( pStar0 * w_phi ) );
 
  double AA = m2_phi - m2;
 
  // construct BW for phi -> rho0 pi0
  double F0 = getF1( m2, m, s12, m2pi0 );
  // double F0 = getF1(m2, m, m2_rho, m2pi0);
  double width0 = getWidth1( m2, m, s12, m2pi0, w_phi );
  // double width0 = getWidth1(m2, m, m2_rho, m2pi0, w_phi);
 
  EvtComplex C0( mw_phi * F0, 0.0 );
  double     BB0  = -m_phi * width0;
  EvtComplex amp0 = C0 / EvtComplex( AA, BB0 );
 
  // construct BW for phi -> rho+ pi-
  double Fp = getF1( m2, m, s13, m2pi );
  // double Fp = getF1(m2, m, m2_rho, m2pi);
  double widthp = getWidth1( m2, m, s13, m2pi, w_phi );
  // double widthp = getWidth1(m2, m, m2_rho, m2pi, w_phi);
 
  EvtComplex Cp( mw_phi * Fp, 0.0 );
  double     BBp  = -m_phi * widthp;
  EvtComplex ampp = Cp / EvtComplex( AA, BBp );
 
  // construct BW for phi -> rho- pi+
  double Fm = getF1( m2, m, s23, m2pi );
  // double Fm = getF1(m2, m, m2_rho, m2pi);
  double widthm = getWidth1( m2, m, s23, m2pi, w_phi );
  // double widthm = getWidth1(m2, m, m2_rho, m2pi, w_phi);
 
  EvtComplex Cm( mw_phi * Fm, 0.0 );
  double     BBm  = -m_phi * widthm;
  EvtComplex ampm = Cm / EvtComplex( AA, BBm );
 
  EvtComplex amp = amp0 + ampp + ampm; // Eq. (15)
 
  double alpham2 = alpha * 2.0;
  F11            = amp * alpham2 * q * H0 * root2;
  F21            = amp * alpham2 * q * ( Hp + Hm );
  F31            = amp * alpham2 * q * ( Hp - Hm );
}
 
double EvtDsTophienu::getF1( double sa, double ma, double sb, double sc ) // a -> b + c
{
  double pStar  = getPStar( sa, sb, sc );
  double pStar0 = getPStar( m2_phi, sb, sc );
  double B      = 1. / sqrt( 1. + rBW2 * pStar * pStar );
  double B0     = 1. / sqrt( 1. + rBW2 * pStar0 * pStar0 );
  double F      = pStar / pStar0 * B / B0;
  return F;
}
 
double EvtDsTophienu::getWidth1( double sa, double ma, double sb, double sc, double width0 ) {
  double pStar  = getPStar( sa, sb, sc );
  double pStar0 = getPStar( m2_phi, sb, sc );
  double F      = getF1( sa, ma, sb, sc );
  double width  = width0 * pStar / pStar0 * m_phi / ma * F * F;
  return width;
}
 
double EvtDsTophienu::getPStar( double s, double s1, double s2 ) {
  double x = s + s1 - s2;
  double t = 0.25 * x * x / s - s1;
  double p;
  if ( t > 0.0 ) { p = sqrt( t ); }
  else
  {
    // std::cout << " Hello, pstar is less than 0.0" << std::endl;
    p = sqrt( -t );
    // p = 0.02;
  }
  return p;
}