EvtD0ToKSKK.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtD0ToKSKK.cc
//
// Description: Routine to handle three-body decays of D0/D0_bar or D+/D- (BESIII arXiv:2006.02800)
//
// Modification history:
//
//    Liaoyuan Dong  Dec 11 05:11:33 2022  Module created
//
//------------------------------------------------------------------------
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include <cmath>
#include <complex>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtD0ToKSKK.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtDecayTable.hh"
using std::endl;
 
EvtD0ToKSKK::~EvtD0ToKSKK() {}
 
void EvtD0ToKSKK::getName(std::string& model_name){
   model_name="D0ToKSKK";
}
 
EvtDecayBase* EvtD0ToKSKK::clone(){
   return new EvtD0ToKSKK;
}
 
void EvtD0ToKSKK::init(){
   checkNArg(0);
   checkNDaug(3);
   checkSpinParent(EvtSpinType::SCALAR);
   checkSpinDaughter(0,EvtSpinType::SCALAR);
   checkSpinDaughter(1,EvtSpinType::SCALAR);
   checkSpinDaughter(2,EvtSpinType::SCALAR);
 
   //cout << "Initializing D0ToKSKK" << endl;
 
   MD = 1.86484; // 0 mass D0
   m1 = 0.4976;  // 1 mass KS
   m2 = 0.49368; // 2 mass K-
   m3 = 0.49368; // 3 mass K+
 
   mag_a00    = 1.0;
   phase_a00  = 0.0;
   spin_a00   = 0;
   mass_a00   = 0.994;   // 23 -> 0:m23sq 8:cosTheta23_12
   gA_a00     = 3.80179;
   gB_a00     = 2.66;
   gC_a00     = 3.80179;
   massB1_a00 = 0.547862;
   massB2_a00 = 0.134977;
   massC1_a00 = 0.497614;
   massC2_a00 = 0.497614;
 
   mag_a0p    = 0.591716;
   phase_a0p  = 2.95792;
   spin_a0p   = 0;
   mass_a0p   = 0.994;  // 13 -> 1:m13sq 5:cosTheta13_12 
   gA_a0p     = 3.80179;
   gB_a0p     = 2.66;
   massB1_a0p = 0.547862;
   massB2_a0p = 0.13957;   
 
   mag_phi    = 0.713488;
   phase_phi  = 1.66557;
   spin_phi   = 1;
   mass_phi   = 1.01946;  // 23 -> 0 8
   width_phi  = 0.004266;
   
   mag_a2p    =  0.126298;
   phase_a2p  = -2.9734;
   spin_a2p   =  2;
   mass_a2p   =  1.3181;  // 13 -> 1 5
   width_a2p  =  0.1098;
   
   mag_a2m    =  0.0965778;
   phase_a2m  = -0.138591;
   spin_a2m   =  2;
   mass_a2m   =  1.3181; // 12 -> 2:m12sq 4:cosTheta12_23
   width_a2m  =  0.1098;
 
   mag_a0_1450m   =  0.210031;
   phase_a0_1450m = -0.23412;
   spin_a0_1450m  =  0;
   mass_a0_1450m  =  1.474; // 12 -> 2 4
   width_a0_1450m =  0.265;
 
   mesonRadius  = 1.5;
   motherRadius = 1.0;
}
 
void EvtD0ToKSKK::initProbMax() {
   setProbMax(143.0);
}
 
void EvtD0ToKSKK::decay( EvtParticle *p){
/*
   double maxprob = 0.0;
   for(int ir=0;ir<=60000000;ir++){
      p->initializePhaseSpace(getNDaug(),getDaugs());
      EvtVector4R D1 = p->getDaug(0)->getP4();
      EvtVector4R D2 = p->getDaug(1)->getP4();
      EvtVector4R D3 = p->getDaug(2)->getP4();
      double value1 = EvaluateAmp(D1, D2, D3);
      if(value1>maxprob) {
         maxprob=value1;
         cout << "ir = " << ir << " maxprob = " << value1 << endl;
      }
   }
   cout << "MAXprob = " << maxprob << endl;
*/
   p->initializePhaseSpace(getNDaug(),getDaugs());
   EvtVector4R D1 = p->getDaug(0)->getP4();  // KS0
   EvtVector4R D2 = p->getDaug(1)->getP4();  // K-
   EvtVector4R D3 = p->getDaug(2)->getP4();  // K+
   double value = EvaluateAmp(D1, D2, D3);
   setProb(value);
   return;
}
 
double EvtD0ToKSKK::EvaluateAmp(EvtVector4R& p4_d1, EvtVector4R& p4_d2, EvtVector4R& p4_d3) {
   double m23sq = (p4_d2 + p4_d3).mass2();
   double m13sq = (p4_d1 + p4_d3).mass2();
   double m12sq = (p4_d1 + p4_d2).mass2();
   double cosTheta23_1v2 = helicityAngle(MD,m2,m3,m1,m23sq,m12sq);  // Angle in RF12 between 3 and 2
   double cosTheta13_1v2 = helicityAngle(MD,m1,m3,m2,m13sq,m12sq);  // Angle in RF13 between 2 and 1
   double cosTheta12_2v3 = helicityAngle(MD,m2,m1,m3,m12sq,m23sq);  // Angle in RF23 between 1 and 2
   //res[0] = m23sq; res[1] = m13sq; res[2] = m12sq; res[3] = cosTheta23_1v2; res[4] = cosTheta13_1v2; res[5] = cosTheta12_2v3;
   
   std::complex<double> Amp_a00 = GetCoefficient(mag_a00,phase_a00)*AmpFlatteRes(m23sq,mass_a00,m2,m3,gA_a00,massB1_a00,massB2_a00,gB_a00,massC1_a00,massC2_a00,gC_a00,spin_a00,mesonRadius)/10.194701039;
   std::complex<double> Amp_a0p = GetCoefficient(mag_a0p,phase_a0p)*AmpFlatteRes(m13sq,mass_a0p,m1,m3,gA_a0p,massB1_a0p,massB2_a0p,gB_a0p,spin_a0p,mesonRadius)/12.861154457;
   std::complex<double> Amp_phi = GetCoefficient(mag_phi,phase_phi)*AmpRelBreitWignerRes(m23sq,mass_phi,m2,m3,width_phi,spin_phi,mesonRadius) *Wigner_d(spin_phi,0,0,acos(cosTheta23_1v2))/715.1406308838;
   std::complex<double> Amp_a2p = GetCoefficient(mag_a2p,phase_a2p)*AmpRelBreitWignerRes(m13sq,mass_a2p,m1,m3,width_a2p,spin_a2p,mesonRadius) *Wigner_d(spin_a2p,0,0,acos(cosTheta13_1v2))/340.70301058;
   std::complex<double> Amp_a2m = GetCoefficient(mag_a2m,phase_a2m)*AmpRelBreitWignerRes(m12sq,mass_a2m,m1,m2,width_a2m,spin_a2m,mesonRadius) *Wigner_d(spin_a2m,0,0,acos(cosTheta12_2v3))/340.737568318;
   std::complex<double> Amp_a0_1450m = GetCoefficient(mag_a0_1450m,phase_a0_1450m)*AmpRelBreitWignerRes(m12sq,mass_a0_1450m,m1,m2,width_a0_1450m,spin_a0_1450m,mesonRadius)/8.9320553176;
   //res[6]  = std::norm(Amp_a00); res[7]  = std::norm(Amp_a0p); res[8]  = std::norm(Amp_phi); res[9]  = std::norm(Amp_a2p); res[10] = std::norm(Amp_a2m); res[11] = std::norm(Amp_a0_1450m);
 
   std::complex<double> Amp = Amp_a00 + Amp_a0p + Amp_phi + Amp_a2p + Amp_a2m + Amp_a0_1450m;
   double value = std::norm(Amp); //res[12] = std::norm(Amp);
   return value;
}
 
double EvtD0ToKSKK::helicityAngle(double M, double m, double m2, double mSpec, double invMassSqA, double invMassSqB) {
   // Calculate energy and momentum of m1/m2 in the invMassSqA rest frame
   double eCms = (invMassSqA + m*m - m2*m2)/(2*sqrt(invMassSqA));
   double pCms = eCms*eCms - m*m;
   // Calculate energy and momentum of mSpec in the invMassSqA rest frame    
   double eSpecCms = (M*M - mSpec*mSpec - invMassSqA)/( 2*sqrt(invMassSqA));
   double pSpecCms = eSpecCms*eSpecCms - mSpec*mSpec;
   double cosAngle = -(invMassSqB - m*m - mSpec*mSpec - 2*eCms*eSpecCms)/(2*sqrt(pCms*pSpecCms));
   return cosAngle;
}
 
std::complex<double> EvtD0ToKSKK::AmpRelBreitWignerRes(double mSq, double mR, double ma, double mb, double width, unsigned int J, double mesonRadius) {
   //J Angular momentum between A&B. In case A&B have spin 0 this is the spin for the resonace.
   //mSq Invariant mass, mR Resonance mass, ma Mass particle A, mb Mass particle B, width  Width of resonance, mesonRadius Scale of interaction range.
   std::complex<double> i(0,1);
   double sqrtS = sqrt(mSq);
   double barrier = FormFactor(sqrtS,ma,mb,J,mesonRadius)/FormFactor(mR,ma,mb,J,mesonRadius);
   std::complex<double> qTerm   = pow( (phspFactor(sqrtS,ma,mb)/phspFactor(mR,ma,mb))*sqrtS/mR, (2*J+1) )*mR/sqrtS;
   std::complex<double> g_final = widthToCoupling(mSq,mR,width,ma,mb,J,mesonRadius);
   std::complex<double> denom   = std::complex<double>(mR*mR-mSq, 0) + (-1.0)*i*sqrtS*(width*qTerm*barrier);
   std::complex<double> result  = g_final/denom;
   return result;
}
 
std::complex<double> EvtD0ToKSKK::widthToCoupling(double mSq, double mR, double width, double ma, double mb, double spin, double mesonRadius) {
   double sqrtS = sqrt(mSq);
   std::complex<double> gammaA(1,0); //spin==0
   if(spin>0){
      std::complex<double> q = qValue(mR,ma,mb);
      double ffR = FormFactor(mR,ma,mb,spin,mesonRadius);
      std::complex<double> qR = pow(q,spin);
      gammaA = ffR*qR;
   }
   std::complex<double> rho = phspFactor(sqrtS,ma,mb);
   std::complex<double> denom = gammaA*sqrt(rho);
   std::complex<double> res = std::complex<double>(sqrt(mR*width), 0)/denom;
   return res;
}
 
std::complex<double> EvtD0ToKSKK::AmpFlatteRes(double mSq, double mR, double massA1, double massA2, double gA, double massB1, double massB2, double couplingB, unsigned int J, double mesonRadius) {
   std::complex<double> i(0,1);
   double sqrtS = sqrt(mSq);
   std::complex<double> termA = gA*gA*phspFactor(sqrtS,massA1,massA2);                 // channel A - signal channel
   std::complex<double> termB = couplingB*couplingB*phspFactor(sqrtS, massB1, massB2); // channel B - hidden channel
   std::complex<double> denom  = std::complex<double>(mR*mR-mSq, 0) + (-1.0)*i*(termA + termB);
   std::complex<double> result = std::complex<double>(gA, 0) / denom;
   return result;
}
std::complex<double> EvtD0ToKSKK::AmpFlatteRes(double mSq, double mR, double massA1, double massA2, double gA, double massB1, double massB2, double couplingB, double massC1, double massC2, double couplingC, unsigned int J, double mesonRadius) {
   std::complex<double> i(0,1);
   double sqrtS = sqrt(mSq);
   std::complex<double> termA = gA*gA*phspFactor(sqrtS,massA1,massA2);                 // channel A - signal channel
   std::complex<double> termB = couplingB*couplingB*phspFactor(sqrtS, massB1, massB2); // channel B - hidden channel
   std::complex<double> termC = couplingC*couplingC*phspFactor(sqrtS, massC1, massC2); // channel C - hidden channel
   std::complex<double> denom  = std::complex<double>(mR*mR-mSq, 0) + (-1.0)*i*(termA + termB + termC);
   std::complex<double> result = std::complex<double>(gA, 0) / denom;
   return result;
}
 
double EvtD0ToKSKK::FormFactor(double sqrtS, double ma, double mb, double spin, double mesonRadius) {
   if(spin==0){ return 1.0; }
   std::complex<double> q = qValue(sqrtS,ma,mb);
   //Blatt-Weisskopt form factors with normalization F(x=mR) = 1. Reference: S.U.Chung Annalen der Physik 4(1995) 404-430
   //z = q / (interaction range). For the interaction range we assume 1/mesonRadius
   double qSq = std::norm(q);
   double z = qSq*mesonRadius*mesonRadius;
   z = fabs(z);
   if (spin==1){
      return( sqrt(2*z/(z+1)) );
   }
   else if (spin==2) {
      return ( sqrt( 13*z*z/( (z-3)*(z-3)+9*z ) ) );
   }
   return 0;
}
 
std::complex<double> EvtD0ToKSKK::phspFactor(double sqrtS, double ma, double mb) {
   // Two types of analytic continuation
   // 1) Complex sqrt
   //      std::complex<double> rhoOld;
   //      rhoOld = sqrt(std::complex<double>(qSqValue(sqrtS,ma,mb))) / (8*M_PI*sqrtS); //PDG definition
   //      rhoOld = sqrt(std::complex<double>(qSqValue(sqrtS,ma,mb))) / (0.5*sqrtS); //BaBar definition
   //      return rhoOld; //complex sqrt
 
   // 2) Correct analytic continuation
   // proper analytic continuation (PDG 2014 - Resonances (47.2.2))
   // I'm not sure of this is correct for the case of two different masses ma and mb.
   // Furthermore we divide by the factor 16*Pi*Sqrt[s]). This is more or less arbitrary
   // and not mentioned in the reference, but it leads to a good agreement between both approaches.
   double s = sqrtS*sqrtS;
   std::complex<double> i(0,1);
   std::complex<double> rho;
   double q = sqrt( fabs(qSqValue(sqrtS,ma,mb)*4/s) );
   if(s<=0){ //below 0
           rho = -q/M_PI*log(fabs((1+q)/(1-q)));
   } else if( 0<s && s<=(ma+mb)*(ma+mb) ){ //below threshold
           rho = ( -2.0*q/M_PI * atan(1/q) )/(i*16.0*M_PI*sqrtS);
   } else if( (ma+mb)*(ma+mb)<s ){//above threshold
           rho = ( -q/M_PI*log( fabs((1+q)/(1-q)) )+i*q )/(i*16.0*M_PI*sqrtS);
   } 
   return rho;
}
 
std::complex<double> EvtD0ToKSKK::qValue(double sqrtS, double ma, double mb) {
   return phspFactor(sqrtS,ma,mb)*8.0*M_PI*sqrtS;
}
 
double EvtD0ToKSKK::qSqValue(double sqrtS, double ma, double mb) {
   double mapb = ma + mb;
   double mamb = ma - mb;
   double xSq = sqrtS*sqrtS;
   double t1 = xSq - mapb*mapb;
   double t2 = xSq - mamb*mamb;
   return ( t1*t2/(4*xSq) );
}
 
std::complex<double> EvtD0ToKSKK::GetCoefficient(double mag, double phase) const {
   return std::complex<double>(fabs(mag)*cos(phase), fabs(mag)*sin(phase));
}
 
double EvtD0ToKSKK::Wigner_d(unsigned int __j, unsigned int __m, unsigned int __n, double __beta) {
   int J = (int)(2.*__j);
   int M = (int)(2.*__m);
   int N = (int)(2.*__n);
   int temp_M, k, k_low, k_hi;
   double const_term = 0.0, sum_term = 0.0, d = 1.0;
   int m_p_n, j_p_m, j_p_n, j_m_m, j_m_n;
   int kmn1, kmn2, jmnk, jmk, jnk;
   double kk;
 
   if (J < 0 || abs (M) > J || abs (N) > J) {
      cout << endl;
      cout << "d: you have entered an illegal number for J, M, N." << endl;
      cout << "Must follow these rules: J >= 0, abs(M) <= J, and abs(N) <= J." 
       << endl;
      cout << "J = " << J <<  " M = " << M <<  " N = " << N << endl;
      return 0.;
   }
  
   if (__beta < 0) {
      __beta = fabs (__beta);
      temp_M = M;
      M = N;
      N = temp_M;
   }
 
   m_p_n = (M + N)/2;
   j_p_m = (J + M)/2;
   j_m_m = (J - M)/2;
   j_p_n = (J + N)/2;
   j_m_n = (J - N)/2;
  
   kk = (double)factorial(j_p_m)*(double)factorial(j_m_m)*(double)factorial(j_p_n)*(double)factorial(j_m_n);
   const_term = pow((-1.0),(j_p_m)) * sqrt(kk);
  
   k_low = ((m_p_n>0) ? m_p_n:0);
   k_hi  = ((j_p_m<j_p_n) ? j_p_m:j_p_n);
 
   for (k = k_low; k <= k_hi; k++) {
      kmn1 = 2*k - (M + N)/2;
      jmnk = J + (M + N)/2 - 2*k;
      jmk  = (J + M)/2 - k;
      jnk  = (J + N)/2 - k;
      kmn2 = k - (M + N)/2;
      sum_term += pow((-1.0), (k)) *((pow(cos(__beta/2.0), kmn1)) *(pow(sin(__beta/2.0), jmnk)))/(factorial(k) *factorial(jmk) *factorial(jnk) *factorial(kmn2));
   }
 
   d = const_term*sum_term*sqrt(2.0*__j+1.0);
   return d;
}
 
int EvtD0ToKSKK::factorial(int i) {
   int f = 1;
   if((i == 0)||(i == 1)) f = 1;
   else{
      while(i > 0){
         f = f*i;
         i--;
      }
   }
   return f;
}