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 "EvtGenBase/EvtPatches.hh"
#include <fstream>
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtResonance.hh"
#include "EvtGenBase/EvtResonance2.hh"
#include "EvtGenModels/EvtDsToEtapipi0.hh"
#include <string>
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtDecayTable.hh"
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;
}