Page MenuHomeHEPForge
Files F19251017

EvtD0gammaDalitz.cpp
No OneTemporary
Actions

Size
12 KB
Referenced Files
None
Subscribers
None

EvtD0gammaDalitz.cpp
View Options

/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
#include "EvtGenModels/EvtD0gammaDalitz.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtDecayTable.hh"
#include "EvtGenBase/EvtFlatte.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 <cmath>
#include <cstdlib>
#include <string>
// Initialize the static variables.
const EvtSpinType::spintype& EvtD0gammaDalitz::_SCALAR = EvtSpinType::SCALAR;
const EvtSpinType::spintype& EvtD0gammaDalitz::_VECTOR = EvtSpinType::VECTOR;
const EvtSpinType::spintype& EvtD0gammaDalitz::_TENSOR = EvtSpinType::TENSOR;
const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_EtaPic =
EvtDalitzReso::EtaPic;
const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_PicPicKK =
EvtDalitzReso::PicPicKK;
const EvtDalitzReso::NumType& EvtD0gammaDalitz::_RBW =
EvtDalitzReso::RBW_CLEO_ZEMACH;
const EvtDalitzReso::NumType& EvtD0gammaDalitz::_GS = EvtDalitzReso::GS_CLEO_ZEMACH;
const EvtDalitzReso::NumType& EvtD0gammaDalitz::_KMAT = EvtDalitzReso::K_MATRIX;
const EvtCyclic3::Pair& EvtD0gammaDalitz::_AB = EvtCyclic3::AB;
const EvtCyclic3::Pair& EvtD0gammaDalitz::_AC = EvtCyclic3::AC;
const EvtCyclic3::Pair& EvtD0gammaDalitz::_BC = EvtCyclic3::BC;
std::string EvtD0gammaDalitz::getName()
{
return "D0GAMMADALITZ";
}
EvtDecayBase* EvtD0gammaDalitz::clone()
{
return new EvtD0gammaDalitz;
}
void EvtD0gammaDalitz::init()
{
// check that there are 0 arguments
checkNArg( 0 );
// Check that this model is valid for the specified decay.
checkNDaug( 3 );
checkSpinParent( _SCALAR );
checkSpinDaughter( 0, _SCALAR );
checkSpinDaughter( 1, _SCALAR );
checkSpinDaughter( 2, _SCALAR );
// Get the values of the EvtId objects from the data files.
readPDGValues();
// Get the EvtId of the D0 and its 3 daughters.
getParentId();
EvtId dau[3];
for ( int index = 0; index < 3; index++ ) {
dau[index] = getDaug( index );
}
// Look for K0bar h+ h-. The order will be K[0SL] h+ h-
for ( int index = 0; index < 3; index++ ) {
if ( ( dau[index] == _K0B ) || ( dau[index] == _KS ) ||
( dau[index] == _KL ) ) {
_d1 = index;
} else if ( ( dau[index] == _PIP ) || ( dau[index] == _KP ) ) {
_d2 = index;
} else if ( ( dau[index] == _PIM ) || ( dau[index] == _KM ) ) {
_d3 = index;
} else {
reportInvalidAndExit();
}
}
// Check if we're dealing with Ks pi pi or with Ks K K.
_isKsPiPi = false;
if ( dau[_d2] == _PIP || dau[_d2] == _PIM ) {
_isKsPiPi = true;
}
}
void EvtD0gammaDalitz::initProbMax()
{
setProbMax( 5200. );
}
void EvtD0gammaDalitz::decay( EvtParticle* part )
{
// Check if the D is from a B+- -> D0 K+- decay with the appropriate model.
EvtParticle* parent =
part->getParent(); // If there are no mistakes, should be B+ or B-.
if ( parent != 0 &&
EvtDecayTable::getInstance()->getDecayFunc( parent )->getName() ==
"BTODDALITZCPK" ) {
EvtId parId = parent->getId();
if ( ( parId == _BP ) || ( parId == _BM ) || ( parId == _B0 ) ||
( parId == _B0B ) ) {
_bFlavor = parId;
} else {
reportInvalidAndExit();
}
} else {
reportInvalidAndExit();
}
// Read the D decay parameters from the B decay model.
// Gamma angle in rad.
double gamma = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg(
0 );
// Strong phase in rad.
double delta = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg(
1 );
// Ratio between B->D0K and B->D0barK
double rB = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 2 );
// Same structure for all of these decays.
part->initializePhaseSpace( getNDaug(), getDaugs() );
EvtVector4R pA = part->getDaug( _d1 )->getP4();
EvtVector4R pB = part->getDaug( _d2 )->getP4();
EvtVector4R pC = part->getDaug( _d3 )->getP4();
// Squared invariant masses.
double mSqAB = ( pA + pB ).mass2();
double mSqAC = ( pA + pC ).mass2();
double mSqBC = ( pB + pC ).mass2();
EvtComplex amp( 1.0, 0.0 );
// Direct and conjugated amplitudes.
EvtComplex ampDir;
EvtComplex ampCnj;
if ( _isKsPiPi ) {
// Direct and conjugated Dalitz points.
EvtDalitzPoint pointDir( _mKs, _mPi, _mPi, mSqAB, mSqBC, mSqAC );
EvtDalitzPoint pointCnj( _mKs, _mPi, _mPi, mSqAC, mSqBC, mSqAB );
// Direct and conjugated amplitudes.
ampDir = dalitzKsPiPi( pointDir );
ampCnj = dalitzKsPiPi( pointCnj );
} else {
// Direct and conjugated Dalitz points.
EvtDalitzPoint pointDir( _mKs, _mK, _mK, mSqAB, mSqBC, mSqAC );
EvtDalitzPoint pointCnj( _mKs, _mK, _mK, mSqAC, mSqBC, mSqAB );
// Direct and conjugated amplitudes.
ampDir = dalitzKsKK( pointDir );
ampCnj = dalitzKsKK( pointCnj );
}
if ( _bFlavor == _BP || _bFlavor == _B0 ) {
amp = ampCnj + rB * exp( EvtComplex( 0., delta + gamma ) ) * ampDir;
} else {
amp = ampDir + rB * exp( EvtComplex( 0., delta - gamma ) ) * ampCnj;
}
vertex( amp );
return;
}
EvtComplex EvtD0gammaDalitz::dalitzKsPiPi( const EvtDalitzPoint& point ) const
{
static const EvtDalitzPlot plot( _mKs, _mPi, _mPi, _mD0 );
EvtComplex amp = 0.;
// This corresponds to relativistic Breit-Wigner distributions. Not K-matrix.
// Defining resonances.
static EvtDalitzReso KStarm( plot, _BC, _AC, _VECTOR, 0.893606, 0.0463407,
_RBW );
static EvtDalitzReso KStarp( plot, _BC, _AB, _VECTOR, 0.893606, 0.0463407,
_RBW );
static EvtDalitzReso rho0( plot, _AC, _BC, _VECTOR, 0.7758, 0.1464, _GS );
static EvtDalitzReso omega( plot, _AC, _BC, _VECTOR, 0.78259, 0.00849, _RBW );
static EvtDalitzReso f0_980( plot, _AC, _BC, _SCALAR, 0.975, 0.044, _RBW );
static EvtDalitzReso f0_1370( plot, _AC, _BC, _SCALAR, 1.434, 0.173, _RBW );
static EvtDalitzReso f2_1270( plot, _AC, _BC, _TENSOR, 1.2754, 0.1851, _RBW );
static EvtDalitzReso K0Starm_1430( plot, _BC, _AC, _SCALAR, 1.459, 0.175,
_RBW );
static EvtDalitzReso K0Starp_1430( plot, _BC, _AB, _SCALAR, 1.459, 0.175,
_RBW );
static EvtDalitzReso K2Starm_1430( plot, _BC, _AC, _TENSOR, 1.4256, 0.0985,
_RBW );
static EvtDalitzReso K2Starp_1430( plot, _BC, _AB, _TENSOR, 1.4256, 0.0985,
_RBW );
static EvtDalitzReso sigma( plot, _AC, _BC, _SCALAR, 0.527699, 0.511861,
_RBW );
static EvtDalitzReso sigma2( plot, _AC, _BC, _SCALAR, 1.03327, 0.0987890,
_RBW );
static EvtDalitzReso KStarm_1680( plot, _BC, _AC, _VECTOR, 1.677, 0.205,
_RBW );
// Adding terms to the amplitude with their corresponding amplitude and phase terms.
amp += EvtComplex( .848984, .893618 );
amp += EvtComplex( -1.16356, 1.19933 ) * KStarm.evaluate( point );
amp += EvtComplex( .106051, -.118513 ) * KStarp.evaluate( point );
amp += EvtComplex( 1.0, 0.0 ) * rho0.evaluate( point );
amp += EvtComplex( -.0249569, .0388072 ) * omega.evaluate( point );
amp += EvtComplex( -.423586, -.236099 ) * f0_980.evaluate( point );
amp += EvtComplex( -2.16486, 3.62385 ) * f0_1370.evaluate( point );
amp += EvtComplex( .217748, -.133327 ) * f2_1270.evaluate( point );
amp += EvtComplex( 1.62128, 1.06816 ) * K0Starm_1430.evaluate( point );
amp += EvtComplex( .148802, .0897144 ) * K0Starp_1430.evaluate( point );
amp += EvtComplex( 1.15489, -.773363 ) * K2Starm_1430.evaluate( point );
amp += EvtComplex( .140865, -.165378 ) * K2Starp_1430.evaluate( point );
amp += EvtComplex( -1.55556, -.931685 ) * sigma.evaluate( point );
amp += EvtComplex( -.273791, -.0535596 ) * sigma2.evaluate( point );
amp += EvtComplex( -1.69720, .128038 ) * KStarm_1680.evaluate( point );
return amp;
}
EvtComplex EvtD0gammaDalitz::dalitzKsKK( const EvtDalitzPoint& point ) const
{
static const EvtDalitzPlot plot( _mKs, _mK, _mK, _mD0 );
// Defining resonances.
static EvtDalitzReso a00_980( plot, _AC, _BC, _SCALAR, 0.999, _RBW, .550173,
.324, _EtaPic );
static EvtDalitzReso phi( plot, _AC, _BC, _VECTOR, 1.01943, .00459319, _RBW );
static EvtDalitzReso a0p_980( plot, _AC, _AB, _SCALAR, 0.999, _RBW, .550173,
.324, _EtaPic );
static EvtDalitzReso f0_1370( plot, _AC, _BC, _SCALAR, 1.350, .265, _RBW );
static EvtDalitzReso a0m_980( plot, _AB, _AC, _SCALAR, 0.999, _RBW, .550173,
.324, _EtaPic );
static EvtDalitzReso f0_980( plot, _AC, _BC, _SCALAR, 0.965, _RBW, .695,
.165, _PicPicKK );
static EvtDalitzReso f2_1270( plot, _AC, _BC, _TENSOR, 1.2754, .1851, _RBW );
static EvtDalitzReso a00_1450( plot, _AC, _BC, _SCALAR, 1.474, .265, _RBW );
static EvtDalitzReso a0p_1450( plot, _AC, _AB, _SCALAR, 1.474, .265, _RBW );
static EvtDalitzReso a0m_1450( plot, _AB, _AC, _SCALAR, 1.474, .265, _RBW );
// Adding terms to the amplitude with their corresponding amplitude and phase terms.
EvtComplex amp( 0., 0. ); // Phase space amplitude.
amp += EvtComplex( 1.0, 0.0 ) * a00_980.evaluate( point );
amp += EvtComplex( -.126314, .188701 ) * phi.evaluate( point );
amp += EvtComplex( -.561428, .0135338 ) * a0p_980.evaluate( point );
amp += EvtComplex( .035, -.00110488 ) * f0_1370.evaluate( point );
amp += EvtComplex( -.0872735, .0791190 ) * a0m_980.evaluate( point );
amp += EvtComplex( 0., 0. ) * f0_980.evaluate( point );
amp += EvtComplex( .257341, -.0408343 ) * f2_1270.evaluate( point );
amp += EvtComplex( -.0614342, -.649930 ) * a00_1450.evaluate( point );
amp += EvtComplex( -.104629, .830120 ) * a0p_1450.evaluate( point );
amp += EvtComplex( 0., 0. ) * a0m_1450.evaluate( point );
return 2.8 *
amp; // Multiply by 2.8 in order to reuse the same probmax as Ks pi pi.
}
void EvtD0gammaDalitz::readPDGValues()
{
// Define the EvtIds.
_BP = EvtPDL::getId( "B+" );
_BM = EvtPDL::getId( "B-" );
_B0 = EvtPDL::getId( "B0" );
_B0B = EvtPDL::getId( "anti-B0" );
_D0 = EvtPDL::getId( "D0" );
_D0B = EvtPDL::getId( "anti-D0" );
_KM = EvtPDL::getId( "K-" );
_KP = EvtPDL::getId( "K+" );
_K0 = EvtPDL::getId( "K0" );
_K0B = EvtPDL::getId( "anti-K0" );
_KL = EvtPDL::getId( "K_L0" );
_KS = EvtPDL::getId( "K_S0" );
_PIM = EvtPDL::getId( "pi-" );
_PIP = EvtPDL::getId( "pi+" );
// Read the relevant masses.
_mD0 = EvtPDL::getMass( _D0 );
_mKs = EvtPDL::getMass( _KS );
_mPi = EvtPDL::getMass( _PIP );
_mK = EvtPDL::getMass( _KP );
}
void EvtD0gammaDalitz::reportInvalidAndExit() const
{
EvtGenReport( EVTGEN_ERROR, "EvtD0gammaDalitz" )
<< "EvtD0gammaDalitz: Invalid mode." << std::endl;
exit( 1 );
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 5:46 AM (1 d, 2 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6560330
Default Alt Text
EvtD0gammaDalitz.cpp (12 KB)

Event Timeline

HEPForge.org · .