EvtDsTof0enu.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: EvtDsTof0enu.cc
//         the necessary file: EvtDsTof0enu.hh
//
// Description: Ds -> f0(980) e+ nu,
//
// Modification history:
//   Liaoyuan Dong   Thu Sep 29 00:12:20 CST 2022   Module created
//------------------------------------------------------------------------
#include "EvtDsTof0enu.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 <stdlib.h>
 
EvtDsTof0enu::~EvtDsTof0enu() {}
 
void EvtDsTof0enu::getName( std::string& model_name ) { model_name = "DsTof0enu"; }
 
EvtDecayBase* EvtDsTof0enu::clone() { return new EvtDsTof0enu; }
 
void EvtDsTof0enu::init() {
  static EvtId DsP = EvtPDL::getId( "D_s+" );
  static EvtId DsM = EvtPDL::getId( "D_s-" );
  static EvtId PIM = EvtPDL::getId( "pi-" );
  static EvtId PIP = EvtPDL::getId( "pi+" );
  static EvtId PI0 = EvtPDL::getId( "pi0" );
  static EvtId KM  = EvtPDL::getId( "K-" );
  static EvtId KP  = EvtPDL::getId( "K+" );
  static EvtId KL  = EvtPDL::getId( "K_L0" );
  static EvtId KS  = EvtPDL::getId( "K_S0" );
  static EvtId EP  = EvtPDL::getId( "e+" );
  static EvtId EM  = EvtPDL::getId( "e-" );
  static EvtId MUP = EvtPDL::getId( "mu+" );
  static EvtId MUM = EvtPDL::getId( "mu-" );
 
  checkNArg( 0 );
  checkNDaug( 4 );
  checkSpinParent( EvtSpinType::SCALAR );
  checkSpinDaughter( 0, EvtSpinType::SCALAR );
  checkSpinDaughter( 1, EvtSpinType::SCALAR );
 
  // std::cout << "EvtDsTof0enu ==> Initialization !" << std::endl;
 
  EvtId parnum = getParentId();
  EvtId d1     = getDaug( 0 );
  EvtId d2     = getDaug( 1 );
  EvtId d3     = getDaug( 2 );
  if ( d3 == EP || d3 == EM )
  {
    if ( parnum == DsP || parnum == DsM )
    {
      if ( d1 == PI0 )
      {
        ProbMax = 33.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (pi0pi0) e+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
      else if ( d1 == PIP || d2 == PIP )
      {
        ProbMax = 33.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (pi-pi+) e+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
      else if ( d1 == KP || d2 == KP )
      {
        ProbMax = 29.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (K-K+) e+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
      else if ( d1 == KS )
      {
        ProbMax = 18.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (KS0KS0) e+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
      else if ( d1 == KL )
      {
        ProbMax = 18.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (KL0KL0) e+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
    }
  }
  if ( d3 == MUP || d3 == MUM )
  {
    if ( parnum == DsP || parnum == DsM )
    {
      if ( d1 == PI0 )
      {
        ProbMax = 33.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (pi0pi0) mu+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
      else if ( d1 == PIP || d2 == PIP )
      {
        ProbMax = 33.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (pi-pi+) mu+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
      else if ( d1 == KP || d2 == KP )
      {
        ProbMax = 29.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (K-K+) mu+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
      else if ( d1 == KS )
      {
        ProbMax = 18.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (KS0KS0) mu+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
      else if ( d1 == KL )
      {
        ProbMax = 18.0;
        // std::cout << "EvtDsTof0enu ==> Ds -> f0 (KL0KL0) mu+ nu : ProMax = " << ProbMax <<
        // std::endl;
      }
    }
  }
 
  ciR = EvtComplex( 1.0, 0.0 );
  ciM = EvtComplex( 0.0, 1.0 );
 
  mPi  = 0.13957;
  mPi0 = 0.1349766;
  mK0  = 0.497611;
  mKa  = 0.493677;
  mEt  = 0.547853;
 
  ///////////// for Ds to f0 e nu
  double mADs = 2.5;
  m2ADs       = mADs * mADs;
  mf0         = 0.9399;
  m2f0        = mf0 * mf0;
  flatte_g1   = 0.199 * mf0;
  flatte_g2   = 3.01 * 0.199 * mf0;
  mDs         = 1.96834;
  m2Ds        = mDs * mDs;
}
 
void EvtDsTof0enu::initProbMax() {
  // cout << "EvtDsTof0enu: setProbMax = " << ProbMax << endl;
  setProbMax( ProbMax );
}
 
void EvtDsTof0enu::decay( EvtParticle* p ) {
  /*
     double maxprob = 0.0;
     for(int ir=0;ir<=60000000;ir++){
        p->initializePhaseSpace(getNDaug(),getDaugs());
        EvtVector4R _pi1 = p->getDaug(0)->getP4();
        EvtVector4R _pi2 = p->getDaug(1)->getP4();
        EvtVector4R _e   = p->getDaug(2)->getP4();
        EvtVector4R _nu  = p->getDaug(3)->getP4();
 
        int pid = EvtPDL::getStdHep(p->getDaug(2)->getId());
        int charm;
        if(pid == 11) charm = -1;
        else charm = 1;
        double m2, q2, cosV, cosL, chi;
        KinVGen(_pi1, _pi2, _e, _nu, charm, m2, q2, cosV, cosL, chi);
        double _prob = calPDF(m2, q2, cosV, cosL, chi);
        if(_prob>maxprob) {
            maxprob=_prob;
            cout << "1 = " << _pi1 << endl;
            cout << "2 = " << _pi2 << endl;
            cout << "3 = " << _e << endl;
            cout << "4 = " << _nu << endl;
            std::cout << "Max PDF = " << ir << " charm= " << charm << " prob= " << _prob <<
     std::endl;
        }
     }
     std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
  */
  p->initializePhaseSpace( getNDaug(), getDaugs() );
  EvtVector4R pi1 = p->getDaug( 0 )->getP4();
  EvtVector4R pi2 = p->getDaug( 1 )->getP4();
  EvtVector4R e   = p->getDaug( 2 )->getP4();
  EvtVector4R nu  = p->getDaug( 3 )->getP4();
 
  int pid = EvtPDL::getStdHep( p->getDaug( 2 )->getId() );
  int charm;
  if ( pid == 11 || pid == 13 ) charm = -1;
  else charm = 1;
  double m2, q2, cosV, cosL, chi;
  KinVGen( pi1, pi2, e, nu, charm, m2, q2, cosV, cosL, chi );
  double prob = calPDF( m2, q2, cosV, cosL, chi );
  setProb( prob );
  return;
}
 
void EvtDsTof0enu::KinVGen( EvtVector4R vp4_K, EvtVector4R vp4_Pi, EvtVector4R vp4_Lep,
                            EvtVector4R vp4_Nu, int charm, double& m2, double& q2,
                            double& cosV, double& cosL, double& chi ) {
  EvtVector4R vp4_KPi = vp4_K + vp4_Pi;
  EvtVector4R vp4_W   = vp4_Lep + vp4_Nu;
  m2                  = vp4_KPi.mass2();
  q2                  = vp4_W.mass2();
 
  EvtVector4R boost;
  boost.set( vp4_KPi.get( 0 ), -vp4_KPi.get( 1 ), -vp4_KPi.get( 2 ), -vp4_KPi.get( 3 ) );
  EvtVector4R vp4_Kp = boostTo( vp4_K, boost );
  cosV               = vp4_Kp.dot( vp4_KPi ) / ( vp4_Kp.d3mag() * vp4_KPi.d3mag() );
 
  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_KPi / vp4_KPi.d3mag();
  EvtVector4R C = vp4_Kp.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 EvtDsTof0enu::calPDF( double m2, double q2, double cosV, double cosL, double chi ) {
  EvtComplex F10   = Resonancef0Flatte( m2, q2 );
  double     I1    = 0.25 * abs2( F10 );
  double     sinL  = sqrt( 1. - cosL * cosL );
  double     sinL2 = sinL * sinL;
  double     cos2L = 1.0 - 2.0 * sinL2;
  double     I     = I1 - I1 * cos2L;
  return I;
}
 
EvtComplex EvtDsTof0enu::Resonancef0Flatte( double m2, double q2 ) {
  double     pPiPi     = getPStar( m2Ds, m2, q2 );
  EvtComplex rhopiPpiM = getrho( m2, mPi );  // rho_pi+pi-
  EvtComplex rhopi0pi0 = getrho( m2, mPi0 ); // rho_pi0pi0
  EvtComplex rhoKPKM   = getrho( m2, mKa );  // rho_K+K-
  EvtComplex rhoK0K0   = getrho( m2, mK0 );  // rho_K0K0
  EvtComplex rhopipi   = ( 2.0 / 3.0 ) * rhopiPpiM + ( 1.0 / 3.0 ) * rhopi0pi0;
  EvtComplex rhoKK     = 0.5 * rhoKPKM + 0.5 * rhoK0K0;
  EvtComplex amp =
      ciR / ( ciR * ( m2f0 - m2 ) - ciM * ( flatte_g1 * rhopipi + flatte_g2 * rhoKK ) );
  EvtComplex F10 = amp * pPiPi * mDs / ( 1.0 - q2 / m2ADs );
  return F10;
}
 
double EvtDsTof0enu::getPStar( double sa, double sb, double sc ) {
  double x = sa + sb - sc;
  double q = 0.25 * x * x / sa - sb;
  double p;
  if ( q > 0.0 ) { p = sqrt( q ); }
  else
  {
    // std::cout << "EvtDToa0enu: pstar is less than 0.0" << std::endl;
    p = 0.01;
  }
  return p;
}
 
EvtComplex EvtDsTof0enu::getrho( double m2, double mX ) {
  EvtComplex rho( 0.0, 0.0 );
  if ( ( 1.0 - 4.0 * mX * mX / m2 ) > 0 ) rho = ciR * sqrt( 1.0 - 4.0 * mX * mX / m2 );
  else rho = ciM * sqrt( 4.0 * mX * mX / m2 - 1.0 );
  return rho;
}