EvtDTopi0pi0enu.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: EvtDTopi0pi0enu.cc
//         the necessary file: EvtDTopi0pi0enu.hh
//
// Description: D -> f0(500) (->pi0 pi0) e+ nu, 
//
// Modification history:
//
//    Liaoyuan Dong  Fri Feb 17 00:56:58 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/EvtDTopi0pi0enu.hh"
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtDecayTable.hh"
#include <stdlib.h>
 
EvtDTopi0pi0enu::~EvtDTopi0pi0enu() {}
 
void EvtDTopi0pi0enu::getName(std::string& model_name){
   model_name="DTopi0pi0enu";
}
 
EvtDecayBase* EvtDTopi0pi0enu::clone(){
   return new EvtDTopi0pi0enu;
}
 
void EvtDTopi0pi0enu::init(){
   checkNArg(0);
   checkNDaug(4);
   checkSpinParent(EvtSpinType::SCALAR);
   checkSpinDaughter(0,EvtSpinType::SCALAR);
   checkSpinDaughter(1,EvtSpinType::SCALAR);
 
   //cout << "Initializing EvtDTopi0pi0enu: ProbMax = 2.48532"<< endl;
 
   mA = 2.42;
   m0_S  = 0.953;
   Dp_mD = 1.86962;
   Pi = atan2(0.0,-1.0);
   mPi = 0.1349766;
   mKa = 0.493677;
   mEt = 0.547853;
}
 
void EvtDTopi0pi0enu::initProbMax() {
   setProbMax(2.48532);
}
 
void EvtDTopi0pi0enu::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(0)->getId());
      int charm;
      if(pid == -211) 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(0)->getId());
   int charm;
   if(pid == -211) 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 EvtDTopi0pi0enu::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 EvtDTopi0pi0enu::calPDF(double m2, double q2, double cosV, double cosL, double chi) {
    double m = sqrt(m2);
    double q = sqrt(q2); 
    EvtComplex F10(0.0,0.0);
    EvtComplex f10(0.0,0.0);
    ResonanceSBugg(m, q, f10);
    F10 = f10;
    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;
}
 
void EvtDTopi0pi0enu::ResonanceSBugg(double m, double q, EvtComplex& F10)
{
    double pKPi = getPStar(Dp_mD, m, q);
    double m2  = m*m;
    double q2  = q*q;
    double mA2 = mA*mA;
   
    double sA    = 0.41*mPi*mPi;
    double mr    = m0_S;         //0.953;
    double mr2   = m0_S*m0_S;    //0.908209;// 0.953*0.953;
    double alpha = 1.3;
    double g4pi  = 0.011;
 
    EvtComplex ciR(1.0, 0.0);
    EvtComplex ciM(0.0, 1.0);
    EvtComplex Gamma1(0.0, 0.0);
    EvtComplex Gamma2(0.0, 0.0);
    EvtComplex Gamma3(0.0, 0.0);
    EvtComplex Gamma4(0.0, 0.0);
 
    Gamma1 = getrho(m2,mPi)*getG1(m2,mr)*(m2-sA)/(mr2-sA);
    Gamma2 = getrho(m2,mKa)*0.6*getG1(m2,mr)*(m2/mr2)*exp(-alpha*fabs(m2-4.0*mKa*mKa));
    Gamma3 = getrho(m2,mEt)*0.2*getG1(m2,mr)*(m2/mr2)*exp(-alpha*fabs(m2-4.0*mEt*mEt));
    if (m > 4*mPi) Gamma4 = ciR*mr*g4pi*(1.0+exp(7.082-2.845*mr2))/(1.0+exp(7.082-2.845*m2));
 
    double AA = mr2-m2-getG1(m2,mr)*(m2-sA)*getZ(m2, mr2)/(mr2-sA);
    EvtComplex amp = ciR/(ciR*AA-ciM*(Gamma1+Gamma2+Gamma3+Gamma4));
    F10 = amp*pKPi*Dp_mD/(1.0-q2/mA2);
}
 
double EvtDTopi0pi0enu::getPStar(double m, double m1, double m2)
{
    double s = m*m;
    double s1 = m1*m1;
    double s2 = m2*m2;
    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 = 0.04;
    }
    return p;
}
 
inline double EvtDTopi0pi0enu::getRho(double m2, double mX)
{
    double rho = 0.0;
    if ((1.0-4.0*mX*mX/m2)>0)
         rho = sqrt(1.0-4.0*mX*mX/m2);
    else rho = 0.0; 
    return rho;
}
 
inline EvtComplex EvtDTopi0pi0enu::getrho(double m2, double mX)
{
    EvtComplex rho(0.0, 0.0);
    EvtComplex ciR(1.0, 0.0);
    EvtComplex ciM(0.0, 1.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;
}
 
inline double EvtDTopi0pi0enu::getG1(double m2, double Mr)
{
    double b1  = 1.302;
    double b2  = 0.340;
    double A   = 2.426;
    double Mr2 = Mr*Mr;// 0.953*0.953;
    double gg1 = Mr*(b1+b2*m2)*exp(-(m2-Mr2)/A);
    return gg1;
}
 
inline double EvtDTopi0pi0enu::getZ(double m2, double Mr2)
{
    double zz = (getRho(m2,mPi)*log((1.0-getRho(m2,mPi))/(1.0+getRho(m2,mPi)))
                -getRho(Mr2,mPi)*log((1.0-getRho(Mr2,mPi))/(1.0+getRho(Mr2,mPi))))/Pi;
 
    return zz;
}