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