EvtEtap2gpipi.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:  EvtEtap2gpipi.cc
//
// Description: etaprime -> gamma pipi without/with box contribution
// see PRL 120, 242003 (2018) and refs. therein
//
// Modification history:
// Dec 18 05:52:33 2022  Liaoyuan Dong  Module created
//
//------------------------------------------------------------------------
 
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtEtap2gpipi.hh"
#include <string>
using namespace std;  //::endl;
 
EvtEtap2gpipi::~EvtEtap2gpipi() {}
 
void EvtEtap2gpipi::getName(std::string& model_name){
  model_name="Etap2gpipi";
}
 
EvtDecayBase* EvtEtap2gpipi::clone(){
  return new EvtEtap2gpipi;
}
 
void EvtEtap2gpipi::init(){
  checkNArg(1);
  checkNDaug(3);
  checkSpinParent(EvtSpinType::SCALAR);
  EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
  _flag=getArg(0);
 
  if (_flag==1) {
     Mrho  = 7.76565e-01;
     Grho  = 1.50839e-01;
     delta = 1.60865e-03;
     argdelta = 6.71184e-02;
     delta2 = 4.15936e-01 ;
     argdelta2 = 2.00844e-02;
     Eetap = -2.13798e+01;
  } else {
     Mrho  = 7.72929e-01;    
     Grho  = 1.50184e-01;   
     delta = 1.58745e-03;    
     argdelta = 6.25729e+00; 
     delta2 = 2.58505e-01;
     argdelta2 = 3.28230e+00;
     Eetap = 0.000000;
  }
  Momega = 0.78265;
  Gomega = 0.00849;
  Mrhop = 1.465;
  Grhop = 0.400;
 
  PI = 3.14159265;
  mpi=0.13957018;
}
 
void EvtEtap2gpipi::initProbMax(){
  if (_flag==1) {
     //std::cout << "Initializing EvtEtap2gpipi... flag= " << _flag << " MaxProb= " << 0.0786 << std::endl;
     setProbMax(0.0786);
  } else {
     //std::cout << "Initializing EvtEtap2gpipi... flag= " << _flag << " MaxProb= " << 0.2760 << std::endl;
     setProbMax(0.2760);
  }
}
 
void EvtEtap2gpipi::decay( EvtParticle *p ){
/*
  double maxprob = 0.0;
  for(int ir=0;ir<=60000000;ir++){
     p->initializePhaseSpace(getNDaug(),getDaugs());
     _pd[0]=p->getDaug(0)->getP4();
     _pd[1]=p->getDaug(1)->getP4();
     _pd[2]=p->getDaug(2)->getP4();
     double Prob = AMPsq();
     if(Prob>maxprob) {
         maxprob=Prob;
         std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
     }
  }
  std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
*/
  p->initializePhaseSpace(getNDaug(),getDaugs());
  _pd[0]=p->getDaug(0)->getP4();
  _pd[1]=p->getDaug(1)->getP4();
  _pd[2]=p->getDaug(2)->getP4();
  double prob = AMPsq();
  setProb(prob);
  return;
}
 
double EvtEtap2gpipi::AMPsq(){
  EvtVector4R prho=_pd[0]+_pd[1];
  EvtVector4R mprho(prho.get(0),-1*prho.get(1),-1*prho.get(2),-1*prho.get(3));
  double m2=(_pd[0]+_pd[1]).mass2(); 
  double m=(_pd[0]+_pd[1]).mass(); 
 
  EvtVector4R pi_rhocms = boostTo(_pd[0],mprho);
  EvtVector4R gam_rhocms = boostTo(_pd[2],mprho);
  EvtVector4R etap_rhocms = boostTo(prho+_pd[2],mprho);
  double coshel_pi = pi_rhocms.dot(etap_rhocms)/pi_rhocms.d3mag()/etap_rhocms.d3mag();
  double ang_part = 1-coshel_pi*coshel_pi;
 
  double kgamma=0;
  if(m<0.95778) kgamma = (0.95778*0.95778-m*m)/(2*0.95778);
  double qpim=0;
  if(m>2*mpi) qpim = sqrt(0.25*m*m - pow(mpi,2));
  
  double qpimrho = sqrt(0.25*Mrho*Mrho - pow(mpi,2));
  double qpimrhop = sqrt(0.25*Mrhop*Mrhop - pow(mpi,2));
 
  double hm = (2./PI)*(qpim/m)*log((m+2*qpim)/(2*mpi));
 
  double hmrho  = (2./PI)*(qpimrho/Mrho)*log((Mrho+2*qpimrho)/(2.*mpi));
  double hmrhop = (2./PI)*(qpimrhop/Mrhop)*log((Mrhop+2*qpimrhop)/(2.*mpi));
 
  double d=(3./PI)*(mpi*mpi/pow(qpimrho,2))*log((Mrho+2*qpimrho)/(2*mpi))+(Mrho/(2*PI*qpimrho))-((mpi*mpi*Mrho)/(PI*pow(qpimrho,3)));
  double drhop=(3./PI)*(mpi*mpi/pow(qpimrhop,2))*log((Mrhop+2*qpimrhop)/(2*mpi))+(Mrhop/(2*PI*qpimrhop))-((mpi*mpi*Mrhop)/(PI*pow(qpimrhop,3)));
 
  double dhds = hmrho*(pow(8*pow(qpimrho,2),-1)-pow(2*Mrho*Mrho,-1))+pow(2*PI*Mrho*Mrho,-1);
  double dhdsrhop = hmrhop*(pow(8*pow(qpimrhop,2),-1)-pow(2*Mrhop*Mrhop,-1))+pow(2*PI*Mrhop*Mrhop,-1);
 
  double fm2 = Grho*(Mrho*Mrho/pow(qpimrho,3))*(qpim*qpim*(hm-hmrho)+(Mrho*Mrho-m*m)*qpimrho*qpimrho*(dhds));
  double fm2rhop = Grhop*(Mrhop*Mrhop/pow(qpimrhop,3))*(qpim*qpim*(hm-hmrhop)+(Mrhop*Mrhop-m*m)*qpimrhop*qpimrhop*(dhdsrhop));
 
  double Grhom = Grho*pow((qpim/qpimrho),3)*(Mrho/m);
  double Grhopm = Grhop*pow((qpim/qpimrhop),3)*(Mrhop/m);
 
  double coefficient = (1./(48*PI*PI*PI))*pow(kgamma,2)*pow(qpim,2);
  double coe1= pow((2*sqrt(48*PI*pow(Mrho,-4))),2);
  
  EvtComplex denomGrhom(Mrho*Mrho-m*m+fm2,-1*Mrho*Grhom);
  EvtComplex real(Mrho*Mrho*(1+d*Grho/Mrho),0);
  EvtComplex BWrho= real/denomGrhom;
  
  EvtComplex denomBWomega(Momega*Momega-m*m,-Momega*Gomega);
  EvtComplex real1(Momega*Momega,0);
  EvtComplex BWomega = real1/denomBWomega;
  EvtComplex cdelta(delta*cos(argdelta),delta*sin(argdelta));
  EvtComplex part2=BWrho*cdelta*BWomega*(m*m)/(Momega*Momega);
 
  EvtComplex denomGrhopm(Mrhop*Mrhop-m*m+fm2rhop,-1*Mrhop*Grhopm);
  EvtComplex realpm(Mrhop*Mrhop*(1+drhop*Grhop/Mrhop),0);
  EvtComplex BWrhop = realpm/denomGrhopm;
  EvtComplex cdelta2(delta2*cos(argdelta2),delta2*sin(argdelta2));
  EvtComplex denom(delta2*cos(argdelta2)+1.0,delta2*sin(argdelta2));
  EvtComplex part3=cdelta2*BWrhop;
  
  EvtComplex eof(Eetap,0);
  EvtComplex amp=(( BWrho* (cdelta*BWomega*(m*m)/(Momega*Momega)+1.0) + cdelta2*BWrhop)/(cdelta2+1)*sqrt(coe1)+ eof);
 
  double total = coefficient* abs2(amp);
  double amp2  = total*ang_part;
  return amp2;
}