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 "EvtEtap2gpipi.hh"
#include "../EvtGenBase/EvtGenKine.hh"
#include "../EvtGenBase/EvtPDL.hh"
#include "../EvtGenBase/EvtParticle.hh"
#include "../EvtGenBase/EvtPatches.hh"
#include "../EvtGenBase/EvtReport.hh"
#include "../EvtGenBase/EvtVector4R.hh"
#include <stdlib.h>
#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;
}