EvtDsToEtapipi0.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: EvtDsToEtapipi0.cc
//
// Description: Routine to handle three-body decays of Ds+/Ds-
//
// Modification history:
//    Liaoyuan Dong    Sun Jun 02 01:23:25 2024    Module created
//------------------------------------------------------------------------
#include "EvtDsToEtapipi0.hh"
#include "../EvtGenBase/EvtConst.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 "../EvtGenBase/EvtResonance.hh"
#include "../EvtGenBase/EvtResonance2.hh"
#include <fstream>
#include <stdlib.h>
#include <string>
using std::endl;
 
EvtDsToEtapipi0::~EvtDsToEtapipi0() {}
 
void EvtDsToEtapipi0::getName( std::string& model_name ) { model_name = "DsToEtapipi0"; }
 
EvtDecayBase* EvtDsToEtapipi0::clone() { return new EvtDsToEtapipi0; }
 
void EvtDsToEtapipi0::init() {
  // check that there are 0 arguments
  checkNArg( 0 );
  checkNDaug( 3 );
  checkSpinParent( EvtSpinType::SCALAR );
  checkSpinDaughter( 0, EvtSpinType::SCALAR );
  checkSpinDaughter( 1, EvtSpinType::SCALAR );
  checkSpinDaughter( 2, EvtSpinType::SCALAR );
 
  phi[1] = 0.783116;
  phi[2] = 2.8811;
  rho[1] = 1.50567;
  rho[2] = 0.845892;
  rho[0] = 1;
  phi[0] = 0;
 
  // cout << "DsToEtapipi0 :" << endl;
  // for (int i=0; i<3; i++) {
  //    cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
  // }
 
  mD           = 1.86486;
  mDs          = 1.9683;
  rRes         = 3.0;
  rD           = 5.0;
  metap        = 0.95778;
  mkstr        = 0.89594;
  mk0          = 0.497614;
  mass_Kaon    = 0.49368;
  mass_Pion    = 0.13957;
  mass_Pi0     = 0.1349766;
  math_pi      = 3.1415926;
  mrho         = 0.77549;
  Grho         = 0.1491;
  ma0          = 0.980;
  Ga0          = 0.0756;
  meta         = 0.547862;
  int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
  for ( int i = 0; i < 4; i++ )
  {
    for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
  }
}
 
void EvtDsToEtapipi0::initProbMax() { setProbMax( 213.5 ); }
 
void EvtDsToEtapipi0::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 P1[4], P2[4], P3[4];
        P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
        P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
        P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
 
        double value;
        value = calDalEva(P1, P2, P3);
        if(value>maxprob) {
           maxprob=value;
           printf("ir = %d, maxprob = %.10f\n", ir, value);
        }
     }
     printf("MAXprob = %.10f\n",maxprob);
  */
  p->initializePhaseSpace( getNDaug(), getDaugs() );
  EvtVector4R D1 = p->getDaug( 0 )->getP4();
  EvtVector4R D2 = p->getDaug( 1 )->getP4();
  EvtVector4R D3 = p->getDaug( 2 )->getP4();
 
  double P1[4], P2[4], P3[4];
  P1[0] = D1.get( 0 );
  P1[1] = D1.get( 1 );
  P1[2] = D1.get( 2 );
  P1[3] = D1.get( 3 );
  P2[0] = D2.get( 0 );
  P2[1] = D2.get( 1 );
  P2[2] = D2.get( 2 );
  P2[3] = D2.get( 3 );
  P3[0] = D3.get( 0 );
  P3[1] = D3.get( 1 );
  P3[2] = D3.get( 2 );
  P3[3] = D3.get( 3 );
 
  double value;
  value = calDalEva( P1, P2, P3 );
  setProb( value );
  return;
}
 
Double_t EvtDsToEtapipi0::calDalEva( Double_t P1[], Double_t P2[], Double_t P3[] ) {
  // pi pi0 eta
  // flag two numbers X, Y, X = spin, Y = is Resonant.
  TComplex PDF[4], cof, pdf, module;
  PDF[0] = Spin_factor( P1, P2, P3, 1, 1 );
  PDF[1] = Spin_factor( P1, P2, P3, 1, 0 );
  PDF[2] = Spin_factor( P1, P3, P2, 0, 31 ) - Spin_factor( P2, P3, P1, 0, 30 );
 
  pdf = TComplex( 0, 0 );
  for ( Int_t i = 0; i < 4; i++ )
  {
    cof = TComplex( rho[i] * TMath::Cos( phi[i] ), rho[i] * TMath::Sin( phi[i] ) );
    pdf += cof * PDF[i];
  }
  module = pdf * TComplex::Conjugate( pdf );
  Double_t value;
  value = module.Re();
  return ( value <= 0 ) ? 1e-20 : value;
}
 
TComplex EvtDsToEtapipi0::Spin_factor( Double_t P1[], Double_t P2[], Double_t P3[], Int_t spin,
                                       Int_t flag ) {
  // D-> R P3, R->P1 P2
  // flag 0, non-res, 1, RBW, 30, Flatte for a0^0(980), 31 Flatte for a0^+(980)
  Double_t R[4], s[3], sp2, sDs, B[2], snk, sck;
  for ( Int_t i = 0; i < 4; i++ ) { R[i] = P1[i] + P2[i]; }
  s[0] = SCADot( R, R );
  s[1] = SCADot( P1, P1 );
  s[2] = SCADot( P2, P2 );
  sp2  = SCADot( P3, P3 );
  sDs  = mDs * mDs;
  snk  = mk0 * mk0;
  sck  = mass_Kaon * mass_Kaon;
  TComplex amp( 0, 0 ), prop;
  Double_t tmp( 0. );
  //-----------for prop-------------------------
  Double_t mass( 0.99 );
  TComplex ci( 0, 1 );
  TComplex rhokk, rhopieta;
  if ( spin == 0 )
  {
    if ( flag == 0 ) prop = TComplex::One();
    if ( flag == 1 ) prop = propagatorRBW( ma0, Ga0, s[0], s[1], s[2], 3.0, 0 );
    if ( flag == 30 )
    {
      rhokk    = Flatte_rhoab( s[0], sck, sck );
      rhopieta = Flatte_rhoab( s[0], s[1], s[2] );
      prop = 1.0 / ( mass * mass - s[0] - ci * ( 0.341 * rhopieta + 0.341 * 0.892 * rhokk ) );
    }
    if ( flag == 31 )
    {
      rhokk    = Flatte_rhoab( s[0], snk, sck );
      rhopieta = Flatte_rhoab( s[0], s[1], s[2] );
      prop = 1.0 / ( mass * mass - s[0] - ci * ( 0.341 * rhopieta + 0.341 * 0.892 * rhokk ) );
    }
    amp = prop;
  }
  else if ( spin == 1 )
  {
    Double_t m_res;
    if ( flag == 0 )
    {
      prop  = TComplex::One();
      m_res = mrho;
    }
    if ( flag == 1 )
    {
      prop =
          propagatorRBW( mrho, Grho, s[0], s[1], s[2], 3.0, 1 ); // only rho for res. 1^- term
      m_res = mrho;
    }
    tmp = 0;
    Double_t T1[4], t1[4];
    calt1( R, P3, T1 );
    calt1( P1, P2, t1 );
    B[0] = barrierNeo( 1, s[0], s[1], s[2], 3.0, m_res ); // spin1 only rho
    B[1] = barrierNeo( 1, sDs, s[0], sp2, 5.0, 1.9683 );
    for ( Int_t i = 0; i < 4; i++ ) { tmp += T1[i] * t1[i] * G[i][i]; }
    amp = tmp * prop * B[0] * B[1];
  }
  else if ( spin == 2 )
  {
    tmp = 0;
    Double_t T2[4][4], t2[4][4];
    calt2( R, P3, T2 );
    calt2( P1, P2, t2 );
    B[0] = barrierNeo( 2, s[0], s[1], s[2], 3.0, mrho ); // wait for changed
    B[1] = barrierNeo( 2, sDs, s[0], sp2, 5.0, 1.9683 );
    for ( Int_t i = 0; i < 4; i++ )
    {
      for ( Int_t j = 0; j < 4; j++ ) { tmp += T2[i][j] * t2[j][i] * G[j][j] * G[i][i]; }
    }
    if ( flag == 0 ) prop = TComplex::One();
    if ( flag == 1 ) prop = propagatorRBW( mrho, Grho, s[0], s[1], s[2], 3.0, 1 );
    // If D wave res. term is found, please use corresponding mass and width for each
    // resonance.
    amp = tmp * prop * B[0] * B[1];
  }
  else { cout << "Only S, P, D wave allowed" << endl; }
  return amp;
}
 
void EvtDsToEtapipi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
  res[0] = a1[0] * a2[0] - a1[1] * a2[1];
  res[1] = a1[1] * a2[0] + a1[0] * a2[1];
}
void EvtDsToEtapipi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
  res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
  res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
}
 
double EvtDsToEtapipi0::SCADot( double a1[4], double a2[4] ) {
  double _cal = 0.;
  _cal        = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
  return _cal;
}
double EvtDsToEtapipi0::barrierNeo( int l, double sa, double sb, double sc, double r,
                                    double mR ) {
  double mR2 = mR * mR;
  double pAB = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
  double pR  = 0.25 * ( mR2 + sb - sc ) * ( mR2 + sb - sc ) / mR2 - sb;
  double zAB = pAB * r * r;
  double zR  = pR * r * r;
  double F   = 0.0;
  if ( l == 0 ) F = 1;
  if ( l == 1 ) F = sqrt( ( 1 + zR ) / ( 1 + zAB ) );
  if ( l == 2 ) F = sqrt( ( 9 + 3 * zAB + zAB * zAB ) / ( 9 + 3 * zAB + zAB * zAB ) );
  return F;
}
//------------------spin-------------------------------------------
void EvtDsToEtapipi0::calt1( double daug1[4], double daug2[4], double t1[4] ) {
  double p, pq;
  double pa[4], qa[4];
  for ( int i = 0; i < 4; i++ )
  {
    pa[i] = daug1[i] + daug2[i];
    qa[i] = daug1[i] - daug2[i];
  }
  p  = SCADot( pa, pa );
  pq = SCADot( pa, qa );
  for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
}
void EvtDsToEtapipi0::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
  double p, r;
  double pa[4], t1[4];
  calt1( daug1, daug2, t1 );
  r = SCADot( t1, t1 );
  for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
  p = SCADot( pa, pa );
  for ( int i = 0; i < 4; i++ )
  {
    for ( int j = 0; j < 4; j++ )
    { t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * ( G[i][j] - pa[i] * pa[j] / p ); }
  }
}
//-------------------prop--------------------------------------------
double EvtDsToEtapipi0::wid( double mass, double sa, double sb, double sc, double r, int l ) {
  double widm = 0.;
  double q    = 0.;
  double q0   = 0.;
  double sa0  = mass * mass;
  double m    = sqrt( sa );
  q           = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
  if ( q < 0 ) q = 1e-16;
  q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
  if ( q0 < 0 ) q0 = 1e-16;
  //  double F = barrier(l,sa,sb,sc,r);
  double z  = q * r * r;
  double z0 = q0 * r * r;
  double F  = 0;
  if ( l == 0 ) F = 1;
  if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
  if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
  double t = q / q0;
  widm     = pow( t, l + 0.5 ) * mass / m * F * F;
  return widm;
}
 
void EvtDsToEtapipi0::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
  double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
  double b[2];
  if ( q > 0 )
  {
    rho[0] = 2 * sqrt( q / sa );
    rho[1] = 0;
  }
  else if ( q < 0 )
  {
    rho[0] = 0;
    rho[1] = 2 * sqrt( -q / sa );
  }
}
 
TComplex EvtDsToEtapipi0::Flatte_rhoab( double sa, double sb, double sc ) {
  double   q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
  TComplex rho;
  TComplex ci( 0, 1 );
  if ( q > 0 ) { rho = 2.0 * sqrt( q / sa ) * ( TComplex::One() ); }
  if ( q < 0 ) { rho = 2.0 * sqrt( -q / sa ) * ci; }
  return rho;
}
 
TComplex EvtDsToEtapipi0::propagatorRBW( double mass, double width, double sa, double sb,
                                         double sc, double r, int l ) {
  TComplex ci( 0, 1 );
  TComplex prop =
      1.0 / ( mass * mass - sa - ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
  return prop;
}