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 "EvtDTopi0pi0enu.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>
 
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;
}