EvtD0Toa0enu.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: EvtD0Toa0enu.cc
//         the necessary file: EvtD0Toa0enu.hh
//
// Description: D -> a0(980) e+ nu,
//
// Modification history:
//
//    Liaoyuan Dong   Sep 29 00:13:16 2022  Module created
//
//------------------------------------------------------------------------
#include "EvtD0Toa0enu.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>
 
EvtD0Toa0enu::~EvtD0Toa0enu() {}
 
void EvtD0Toa0enu::getName( std::string& model_name ) { model_name = "D0Toa0enu"; }
 
EvtDecayBase* EvtD0Toa0enu::clone() { return new EvtD0Toa0enu; }
 
void EvtD0Toa0enu::init() {
  checkNArg( 0 );
  checkNDaug( 4 );
  checkSpinParent( EvtSpinType::SCALAR );
  checkSpinDaughter( 0, EvtSpinType::SCALAR );
  checkSpinDaughter( 1, EvtSpinType::SCALAR );
 
  mode    = 0; // D0
  ProbMax = 15.55;
 
  // std::cout << "Initializing EvtD0Toa0enu (D0 -> a0(980)- e+ nu_e) ProbMax = " << ProbMax <<
  // std::endl;
 
  double m_Kaon = 0.493677;
  m2_Kaon       = m_Kaon * m_Kaon;
  double m_K0   = 0.497614;
  m2_K0         = m_K0 * m_K0;
  double m_eta  = 0.547853;
  m2_eta        = m_eta * m_eta;
  double m_Pion = 0.13957;
  m2_Pion       = m_Pion * m_Pion;
  double m_Pi0  = 0.1349766;
  m2_Pi0        = m_Pi0 * m_Pi0;
  m_D0          = 1.86486;
  m2_D0         = m_D0 * m_D0;
  m_D           = 1.86962;
  m2_D          = m_D * m_D;
 
  ciR = EvtComplex( 1.0, 0.0 );
  ciM = EvtComplex( 0.0, 1.0 );
 
  ///////////// for D to a0 e nu
  double mAD0 = 2.42;
  m2AD0       = mAD0 * mAD0;
  double mAD  = 2.42;
  m2AD        = mAD * mAD;
  double m_a0 = 0.990;
  m2_a0       = m_a0 * m_a0;
  flatte_g1   = 0.341;
  flatte_g2   = 0.341 * 0.892;
}
 
void EvtD0Toa0enu::initProbMax() { setProbMax( ProbMax ); }
 
void EvtD0Toa0enu::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;
            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 EvtD0Toa0enu::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 EvtD0Toa0enu::calPDF( double m2, double q2, double cosV, double cosL, double chi ) {
  EvtComplex F10( 0.0, 0.0 );
  if ( mode == 0 ) F10 = a0pFlatte( m2, q2 );
  else if ( mode == 1 ) F10 = a0nFlatte( 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 EvtD0Toa0enu::a0nFlatte( double m2, double q2 ) {
  double     PetaPi   = getPStar( m2_D, m2, q2 );
  EvtComplex rhokk    = Flatte_rhoab( m2, m2_Kaon, m2_Kaon );
  EvtComplex rhopieta = Flatte_rhoab( m2, m2_Pi0, m2_eta );
  EvtComplex amp =
      ciR / ( ciR * ( m2_a0 - m2 ) - ciM * ( flatte_g1 * rhopieta + flatte_g2 * rhokk ) );
  EvtComplex F10 = amp * PetaPi * m_D / ( 1.0 - q2 / m2AD );
  return F10;
}
 
EvtComplex EvtD0Toa0enu::a0pFlatte( double m2, double q2 ) {
  double     PetaPi   = getPStar( m2_D0, m2, q2 );
  EvtComplex rhokk    = Flatte_rhoab( m2, m2_K0, m2_Kaon );
  EvtComplex rhopieta = Flatte_rhoab( m2, m2_Pion, m2_eta );
  EvtComplex amp =
      ciR / ( ciR * ( m2_a0 - m2 ) - ciM * ( flatte_g1 * rhopieta + flatte_g2 * rhokk ) );
  EvtComplex F10 = amp * PetaPi * m_D0 / ( 1.0 - q2 / m2AD0 );
  return F10;
}
 
EvtComplex EvtD0Toa0enu::Flatte_rhoab( double sa, double sb, double sc ) {
  EvtComplex rho( 0.0, 0.0 );
  double     x = sa + sb - sc;
  double     q = 0.25 * x * x / sa - sb;
  if ( q > 0 ) { rho = 2.0 * sqrt( q / sa ) * ciR; }
  else { rho = 2.0 * sqrt( -q / sa ) * ciM; }
  return rho;
}
 
double EvtD0Toa0enu::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 << "EvtD0Toa0enu: pstar is less than 0.0" << std::endl;
    p = 0.01;
  }
  return p;
}