Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenModels/EvtbTosllNPR.cpp
- This file was added.
/*********************************************************************** | |||||
* 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/EvtbTosllNPR.hh" | |||||
#include "EvtGenBase/EvtGenKine.hh" | |||||
#include "EvtGenBase/EvtPDL.hh" | |||||
#include "EvtGenBase/EvtParticle.hh" | |||||
#include "EvtGenBase/EvtScalarParticle.hh" | |||||
#include "EvtGenModels/EvtbTosllAmp.hh" | |||||
#include "EvtGenModels/EvtbTosllBSZFF.hh" | |||||
#include "EvtGenModels/EvtbTosllVectorAmpNP.hh" | |||||
#include <algorithm> | |||||
#include <ccomplex> | |||||
using std::endl; | |||||
std::string EvtbTosllNPR::getName() | |||||
{ | |||||
return "BTOSLLNPR"; | |||||
} | |||||
EvtDecayBase* EvtbTosllNPR::clone() | |||||
{ | |||||
return new EvtbTosllNPR; | |||||
} | |||||
void EvtbTosllNPR::decay( EvtParticle* p ) | |||||
{ | |||||
static EvtVector4R p4[3]; | |||||
static double mass[3]; | |||||
double m_b = p->mass(); | |||||
for ( int i = 0; i < 3; i++ ) | |||||
mass[i] = p->getDaug( i )->mass(); | |||||
EvtId* daughters = getDaugs(); | |||||
double weight = EvtGenKine::PhaseSpacePole2( m_b, mass[2], mass[1], mass[0], | |||||
p4, _ls ); | |||||
p->getDaug( 0 )->init( daughters[0], p4[2] ); | |||||
p->getDaug( 1 )->init( daughters[1], p4[1] ); | |||||
p->getDaug( 2 )->init( daughters[2], p4[0] ); | |||||
setWeight( weight ); | |||||
_calcamp->CalcAmp( p, _amp2, _ffmodel.get() ); | |||||
} | |||||
struct ABW_t { | |||||
double Mv2, MGv, A; | |||||
ABW_t( double m, double gtot, double, double bll ) | |||||
{ | |||||
const double alpha = 1 / 137.035999084; // (+-21) PDG 2021 | |||||
Mv2 = m * m; | |||||
MGv = m * gtot; | |||||
A = 3 * M_PI / ( alpha * alpha ) * bll * gtot; | |||||
} | |||||
std::complex<double> R( double s ) const | |||||
{ | |||||
return A / std::complex<double>( s - Mv2, MGv ); | |||||
} | |||||
void addknots( std::vector<double>& en ) const | |||||
{ | |||||
double Mv = sqrt( Mv2 ), Gv = MGv / Mv; | |||||
double w = 0.6; | |||||
double smin = ( Mv - w * Gv ) * ( Mv - w * Gv ); | |||||
double smax = ( Mv + w * Gv ) * ( Mv + w * Gv ); | |||||
en.push_back( Mv * Mv ); | |||||
for ( int i = 0, n = 30; i < n; i++ ) { | |||||
double s = smin + ( smax - smin ) * ( i + 0.5 ) / n; | |||||
en.push_back( s ); | |||||
} | |||||
w = 2; | |||||
double smin1 = ( Mv - w * Gv ) * ( Mv - w * Gv ); | |||||
double smax1 = ( Mv + w * Gv ) * ( Mv + w * Gv ); | |||||
for ( int i = 0, n = 30; i < n; i++ ) { | |||||
double s = smin1 + ( smax1 - smin1 ) * ( i + 0.5 ) / n; | |||||
if ( s >= smin && s <= smax ) | |||||
continue; | |||||
en.push_back( s ); | |||||
} | |||||
w = 8; | |||||
double smin2 = ( Mv - w * Gv ) * ( Mv - w * Gv ); | |||||
double smax2 = ( Mv + w * Gv ) * ( Mv + w * Gv ); | |||||
for ( int i = 0, n = 30; i < n; i++ ) { | |||||
double s = smin2 + ( smax2 - smin2 ) * ( i + 0.5 ) / n; | |||||
if ( s >= smin1 && s <= smax1 ) | |||||
continue; | |||||
en.push_back( s ); | |||||
} | |||||
w = 30; | |||||
double smin3 = ( Mv - w * Gv ) * ( Mv - w * Gv ); | |||||
double smax3 = ( Mv + w * Gv ) * ( Mv + w * Gv ); | |||||
for ( int i = 0, n = 30; i < n; i++ ) { | |||||
double s = smin3 + ( smax3 - smin3 ) * ( i + 0.5 ) / n; | |||||
if ( s >= smin2 && s <= smax2 ) | |||||
continue; | |||||
en.push_back( s ); | |||||
} | |||||
w = 100; | |||||
double smin4 = ( Mv - w * Gv ) * ( Mv - w * Gv ); | |||||
double smax4 = ( Mv + w * Gv ) * ( Mv + w * Gv ); | |||||
for ( int i = 0, n = 20; i < n; i++ ) { | |||||
double s = smin4 + ( smax4 - smin4 ) * ( i + 0.5 ) / n; | |||||
if ( s >= smin3 && s <= smax3 ) | |||||
continue; | |||||
en.push_back( s ); | |||||
} | |||||
} | |||||
}; | |||||
ABW_t ajpsi( 3.0969, 0.0926 / 1000, 0.0926 / 1000 * 0.877, 0.05971 ); | |||||
ABW_t apsi2s( 3.6861, 0.294 / 1000, 0.294 / 1000 * 0.9785, 0.00793 ); | |||||
ABW_t apsi3770( 3773.7 / 1000, 27.2 / 1000, 27.2 / 1000, 1e-5 ); | |||||
ABW_t apsi4040( 4039 / 1000., 80 / 1000., 80 / 1000., 1.07e-5 ); | |||||
ABW_t apsi4160( 4191 / 1000., 70 / 1000., 70 / 1000., 6.9e-6 ); | |||||
ABW_t apsi4230( 4220 / 1000., 60 / 1000., 60 / 1000., 1e-5 ); | |||||
std::complex<double> aR( double s ) | |||||
{ | |||||
return ajpsi.R( s ) + apsi2s.R( s ) + apsi3770.R( s ) + apsi4040.R( s ) + | |||||
apsi4160.R( s ) + apsi4230.R( s ); | |||||
} | |||||
std::vector<double> getgrid() | |||||
{ | |||||
const double m_c = 1.3, m_s = 0.2, m_e = 0.511e-3, MD0 = 1864.84 / 1000; | |||||
std::vector<ABW_t> v = { ajpsi, apsi2s, apsi3770, | |||||
apsi4040, apsi4160, apsi4230 }; | |||||
std::vector<double> en; | |||||
for ( const auto& t : v ) | |||||
t.addknots( en ); | |||||
double smax = 25; //4.37*4.37; | |||||
for ( unsigned i = 0, n = 200; i < n; i++ ) { | |||||
double s = smax / n * ( i + 0.5 ); | |||||
en.push_back( s ); | |||||
} | |||||
en.push_back( ajpsi.Mv2 + 0.1 ); | |||||
en.push_back( ajpsi.Mv2 - 0.1 ); | |||||
en.push_back( ajpsi.Mv2 + 0.08 ); | |||||
en.push_back( ajpsi.Mv2 - 0.08 ); | |||||
en.push_back( ajpsi.Mv2 + 0.065 ); | |||||
en.push_back( ajpsi.Mv2 - 0.065 ); | |||||
double t0; | |||||
for ( t0 = 4 * m_e * m_e; t0 < 0.15; t0 *= 1.5 ) { | |||||
en.push_back( t0 ); | |||||
} | |||||
for ( t0 = 4 * m_e * m_e; t0 < 0.5; t0 *= 1.5 ) { | |||||
en.push_back( 4 * MD0 * MD0 + t0 ); | |||||
en.push_back( 4 * MD0 * MD0 - t0 ); | |||||
} | |||||
en.push_back( 4 * m_c * m_c ); | |||||
for ( t0 = 0.00125; t0 < 0.05; t0 *= 1.5 ) { | |||||
en.push_back( 4 * m_c * m_c + t0 ); | |||||
en.push_back( 4 * m_c * m_c - t0 ); | |||||
} | |||||
en.push_back( 4 * m_s * m_s ); | |||||
for ( t0 = 0.00125; t0 < 0.1; t0 *= 1.5 ) { | |||||
en.push_back( 4 * m_s * m_s + t0 ); | |||||
en.push_back( 4 * m_s * m_s - t0 ); | |||||
} | |||||
std::sort( en.begin(), en.end() ); | |||||
for ( std::vector<double>::iterator it = en.begin(); it != en.end(); it++ ) { | |||||
if ( *it > smax ) { | |||||
en.erase( it, en.end() ); | |||||
break; | |||||
} | |||||
} | |||||
return en; | |||||
} | |||||
void EvtbTosllNPR::initProbMax() | |||||
{ | |||||
EvtId pId = getParentId(), mId = getDaug( 0 ), l1Id = getDaug( 1 ), | |||||
l2Id = getDaug( 2 ); | |||||
// //This routine sets the _poleSize. | |||||
// double mymaxprob = _calcamp->CalcMaxProb( pId, mId, l1Id, l2Id, _ffmodel.get(), _poleSize ); | |||||
// setProbMax( mymaxprob ); | |||||
std::vector<double> v = getgrid(); | |||||
std::vector<EvtPointf> pp; | |||||
EvtScalarParticle* scalar_part; | |||||
EvtParticle* root_part; | |||||
scalar_part = new EvtScalarParticle; | |||||
//cludge to avoid generating random numbers! | |||||
scalar_part->noLifeTime(); | |||||
EvtVector4R p_init; | |||||
p_init.set( EvtPDL::getMass( pId ), 0.0, 0.0, 0.0 ); | |||||
scalar_part->init( pId, p_init ); | |||||
root_part = (EvtParticle*)scalar_part; | |||||
root_part->setDiagonalSpinDensity(); | |||||
EvtParticle *daughter, *lep1, *lep2; | |||||
EvtAmp amp; | |||||
EvtId listdaug[3]; | |||||
listdaug[0] = mId; | |||||
listdaug[1] = l1Id; | |||||
listdaug[2] = l2Id; | |||||
amp.init( pId, 3, listdaug ); | |||||
root_part->makeDaughters( 3, listdaug ); | |||||
daughter = root_part->getDaug( 0 ); | |||||
lep1 = root_part->getDaug( 1 ); | |||||
lep2 = root_part->getDaug( 2 ); | |||||
//cludge to avoid generating random numbers! | |||||
daughter->noLifeTime(); | |||||
lep1->noLifeTime(); | |||||
lep2->noLifeTime(); | |||||
//Initial particle is unpolarized, well it is a scalar so it is trivial | |||||
EvtSpinDensity rho; | |||||
rho.setDiag( root_part->getSpinStates() ); | |||||
double mass[3]; | |||||
double m = root_part->mass(); | |||||
EvtVector4R p4meson, p4lepton1, p4lepton2, p4w; | |||||
double maxprobfound = 0.0; | |||||
mass[1] = EvtPDL::getMeanMass( l1Id ); | |||||
mass[2] = EvtPDL::getMeanMass( l2Id ); | |||||
std::vector<double> mH; | |||||
mH.push_back( EvtPDL::getMeanMass( mId ) ); | |||||
//if the particle is narrow dont bother with changing the mass. | |||||
double g = EvtPDL::getWidth( mId ); | |||||
if ( g > 0 ) { | |||||
mH.push_back( EvtPDL::getMinMass( mId ) ); | |||||
mH.push_back( | |||||
std::min( EvtPDL::getMaxMass( mId ), m - mass[1] - mass[2] ) ); | |||||
mH.push_back( mH.front() - g ); | |||||
mH.push_back( mH.front() + g ); | |||||
} | |||||
double q2min = ( mass[1] + mass[2] ) * ( mass[1] + mass[2] ); | |||||
double m0 = EvtPDL::getMinMass( mId ); | |||||
double q2max = ( m - m0 ) * ( m - m0 ); | |||||
m0 = mH[0]; | |||||
//loop over q2 | |||||
for ( double q2 : v ) { | |||||
// want to avoid picking up the tail of the photon propagator | |||||
if ( !( q2min <= q2 && q2 < q2max ) ) | |||||
continue; | |||||
double Erho = ( m * m + m0 * m0 - q2 ) / ( 2.0 * m ), | |||||
Prho = sqrt( ( Erho - m0 ) * ( Erho + m0 ) ); | |||||
p4meson.set( Erho, 0.0, 0.0, -Prho ); | |||||
p4w.set( m - Erho, 0.0, 0.0, Prho ); | |||||
//This is in the W rest frame | |||||
double Elep = sqrt( q2 ) * 0.5, | |||||
Plep = sqrt( Elep * Elep - mass[1] * mass[1] ); | |||||
const int nj = 3 + 2 + 4 + 8 + 32; | |||||
double cmin = -1, cmax = 1, dc = ( cmax - cmin ) / ( nj - 1 ); | |||||
double maxprob = 0; | |||||
for ( int j = 0; j < nj; j++ ) { | |||||
double c = cmin + dc * j; | |||||
//These are in the W rest frame. Need to boost out into the B frame. | |||||
double Py = Plep * sqrt( 1.0 - c * c ), Pz = Plep * c; | |||||
p4lepton1.set( Elep, 0.0, Py, Pz ); | |||||
p4lepton2.set( Elep, 0.0, -Py, -Pz ); | |||||
p4lepton1 = boostTo( p4lepton1, p4w ); | |||||
p4lepton2 = boostTo( p4lepton2, p4w ); | |||||
//Now initialize the daughters... | |||||
daughter->init( mId, p4meson ); | |||||
lep1->init( l1Id, p4lepton1 ); | |||||
lep2->init( l2Id, p4lepton2 ); | |||||
_calcamp->CalcAmp( root_part, amp, _ffmodel.get() ); | |||||
double prob = rho.normalizedProb( amp.getSpinDensity() ); | |||||
maxprob = std::max( maxprob, prob ); | |||||
} | |||||
// printf("%f %f\n",q2,maxprob); | |||||
pp.push_back( { (float)q2, (float)maxprob } ); | |||||
maxprobfound = std::max( maxprobfound, maxprob ); | |||||
} | |||||
root_part->deleteTree(); | |||||
_ls.init( pp ); | |||||
maxprobfound *= 8e-8 * 1.1; | |||||
setProbMax( maxprobfound ); | |||||
} | |||||
void EvtbTosllNPR::init() | |||||
{ | |||||
// First choose format of NP coefficients from the .DEC file | |||||
// Cartesian(x,y)(0) or Polar(R,phase)(1) | |||||
int n = getNArg(); | |||||
checkNDaug( 3 ); | |||||
//We expect the parent to be a scalar | |||||
//and the daughters to be K* lepton+ lepton- | |||||
checkSpinParent( EvtSpinType::SCALAR ); | |||||
checkSpinDaughter( 1, EvtSpinType::DIRAC ); | |||||
checkSpinDaughter( 2, EvtSpinType::DIRAC ); | |||||
EvtId mId = getDaug( 0 ); | |||||
EvtId IdKst0 = EvtPDL::getId( "K*0" ), IdaKst0 = EvtPDL::getId( "anti-K*0" ), | |||||
IdKstp = EvtPDL::getId( "K*+" ), IdKstm = EvtPDL::getId( "K*-" ); | |||||
EvtId IdK0 = EvtPDL::getId( "K0" ), IdaK0 = EvtPDL::getId( "anti-K0" ), | |||||
IdKp = EvtPDL::getId( "K+" ), IdKm = EvtPDL::getId( "K-" ); | |||||
if ( mId != IdKst0 && mId != IdaKst0 && mId != IdKstp && mId != IdKstm && | |||||
mId != IdK0 && mId != IdaK0 && mId != IdKp && mId != IdKm ) { | |||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | |||||
<< "EvtbTosllNPR generator expected a K(*) 1st daughter, found: " | |||||
<< EvtPDL::name( getDaug( 0 ) ) << endl; | |||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | |||||
<< "Will terminate execution!" << endl; | |||||
::abort(); | |||||
} | |||||
_ffmodel = std::make_unique<EvtbTosllBSZFF>(); | |||||
_calcamp = std::make_unique<EvtbTosllVectorAmpNP>(); | |||||
auto getInteger = [this]( int i ) -> int { | |||||
double a = getArg( i ); | |||||
if ( a - static_cast<int>( a ) != 0 ) { | |||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | |||||
<< "Parameters is not integer in the BTOSLLNPR decay model: " << i | |||||
<< " " << a << endl; | |||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | |||||
<< "Will terminate execution!" << endl; | |||||
::abort(); | |||||
} | |||||
return static_cast<int>( a ); | |||||
}; | |||||
EvtbTosllVectorAmpNP* amp = static_cast<EvtbTosllVectorAmpNP*>( | |||||
_calcamp.get() ); | |||||
double phi = 0; | |||||
if ( n > 0 ) { // parse arguments | |||||
int i = 0, coordsyst = getInteger( i++ ); | |||||
auto getcomplex = [this, &i, coordsyst]() -> EvtComplex { | |||||
double a0 = getArg( i++ ); | |||||
double a1 = getArg( i++ ); | |||||
return ( coordsyst ) ? EvtComplex( a0 * cos( a1 ), a0 * sin( a1 ) ) | |||||
: EvtComplex( a0, a1 ); | |||||
}; | |||||
auto getreal = [this, &i]() -> double { return getArg( i++ ); }; | |||||
while ( i < n ) { | |||||
int id = getInteger( i++ ); // New Physics cooeficient Id | |||||
if ( id == 0 ) | |||||
amp->m_dc7 = getcomplex(); // delta C_7eff | |||||
if ( id == 1 ) | |||||
amp->m_dc9 = getcomplex(); // delta C_9eff | |||||
if ( id == 2 ) | |||||
amp->m_dc10 = getcomplex(); // delta C_10eff | |||||
if ( id == 3 ) | |||||
amp->m_c7p = getcomplex(); // C'_7eff -- right hand polarizations | |||||
if ( id == 4 ) | |||||
amp->m_c9p = getcomplex(); // C'_9eff -- right hand polarizations | |||||
if ( id == 5 ) | |||||
amp->m_c10p = getcomplex(); // c'_10eff -- right hand polarizations | |||||
if ( id == 6 ) | |||||
amp->m_cS = getcomplex(); // (C_S - C'_S) -- scalar right and left polarizations | |||||
if ( id == 7 ) | |||||
amp->m_cP = getcomplex(); // (C_P - C'_P) -- pseudo-scalar right and left polarizations | |||||
if ( id == 10 ) | |||||
phi = getreal(); | |||||
} | |||||
} | |||||
std::vector<double> v = getgrid(); | |||||
std::complex<double> eiphi = std::complex<double>( cos( phi ), sin( phi ) ); | |||||
for ( auto t : v ) | |||||
amp->m_reso.push_back( std::make_pair( t, -aR( t ) * eiphi ) ); | |||||
} |