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 "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtDsTophienu.hh"
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtDecayTable.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;
}