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