EvtEtap2pipieta.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:  EvtEtap2pipieta.cc
//
// Description: etaprime -> pi+pi- eta
// see PHYSICAL REVIEW D 97, 012003 (2018)
//
// Modification history:
// Jun 03 00:42:10 2024  Liaoyuan Dong  Module created
//------------------------------------------------------------------------
 
#include "EvtEtap2pipieta.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;
 
EvtEtap2pipieta::~EvtEtap2pipieta() {}
 
void EvtEtap2pipieta::getName( std::string& model_name ) { model_name = "Etap2pipieta"; }
 
EvtDecayBase* EvtEtap2pipieta::clone() { return new EvtEtap2pipieta; }
 
void EvtEtap2pipieta::init() {
  checkNDaug( 3 );
  checkSpinParent( EvtSpinType::SCALAR );
  EvtSpinType::spintype parenttype = EvtPDL::getSpinType( getParentId() );
  m_MetaOverMpip                   = 0.547862 / 0.13957018; // m_Meta/m_Mpi
  m_sqrt3                          = sqrt( 3.0 );
  par[0]                           = -0.127;
  par[1]                           = -0.106;
  par[2]                           = 1.5e-03;
  par[3]                           = -0.082;
}
 
void EvtEtap2pipieta::initProbMax() { setProbMax( 1.03805 ); }
 
void EvtEtap2pipieta::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 EvtEtap2pipieta::AMPsq() {
  EvtVector4R petap = _pd[0] + _pd[1] + _pd[2];
  EvtVector4R Betap( petap.get( 0 ), -1.0 * petap.get( 1 ), -1.0 * petap.get( 2 ),
                     -1.0 * petap.get( 3 ) );
 
  EvtVector4R pip_etapcms = boostTo( _pd[0], Betap );
  EvtVector4R pim_etapcms = boostTo( _pd[1], Betap );
  EvtVector4R eta_etapcms = boostTo( _pd[2], Betap );
 
  double Q = Betap.mass() - pip_etapcms.mass() - pim_etapcms.mass() - eta_etapcms.mass();
 
  double Tp = pip_etapcms.get( 0 ) - pip_etapcms.mass();
  double Tm = pim_etapcms.get( 0 ) - pim_etapcms.mass();
  double Te = eta_etapcms.get( 0 ) - eta_etapcms.mass();
 
  double X = m_sqrt3 * ( Tp - Tm ) / Q;
  double Y = ( 2.0 + m_MetaOverMpip ) * Te / Q - 1.0;
 
  double amp2 = 1.0 + par[0] * Y + par[1] * Y * Y + par[2] * X + par[3] * X * X;
  return amp2;
}