Page MenuHomeHEPForge
Differential D97 Diff 415 src/EvtGenModels/EvtVtoSll.cpp

Changeset View

Standalone View

View Options

src/EvtGenModels/EvtVtoSll.cpp

Show All 16 Lines
* You should have received a copy of the GNU General Public License * * You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. * * along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/ ***********************************************************************/
#include "EvtGenModels/EvtVtoSll.hh" #include "EvtGenModels/EvtVtoSll.hh"
#include "EvtGenBase/EvtDiracSpinor.hh" #include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtGenKine.hh" #include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtIdSet.hh"
#include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtTensor4C.hh" #include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh" #include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtVectorParticle.hh"
#include <iostream> #include <iostream>
#include <stdlib.h> #include <stdlib.h>
#include <string> #include <string>
std::string EvtVtoSll::getName() std::string EvtVtoSll::getName()
{ {
return "VTOSLL"; return "VTOSLL";
Show All 10 Linesvoid EvtVtoSll::init()
checkNArg( 0 ); checkNArg( 0 );
checkNDaug( 3 ); checkNDaug( 3 );
checkSpinParent( EvtSpinType::VECTOR ); checkSpinParent( EvtSpinType::VECTOR );
checkSpinDaughter( 0, EvtSpinType::SCALAR ); checkSpinDaughter( 0, EvtSpinType::SCALAR );
checkSpinDaughter( 1, EvtSpinType::DIRAC ); checkSpinDaughter( 1, EvtSpinType::DIRAC );
checkSpinDaughter( 2, EvtSpinType::DIRAC ); checkSpinDaughter( 2, EvtSpinType::DIRAC );
// Work out whether we have electron mode
const EvtIdSet leptons{ "e-", "e+" };
if ( leptons.contains( getDaug( 1 ) ) || leptons.contains( getDaug( 2 ) ) ) {
m_electronMode = true;
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< " EvtVtoSll has dielectron final state" << std::endl;
}
} }
void EvtVtoSll::initProbMax() void EvtVtoSll::initProbMax()
{ {
setProbMax( 10.0 ); EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< " EvtVtoSll is finding maximum probability ... " << std::endl;
double theProbMax = 0;
double theProbMax_q2 = 0;
double theProbMax_theta = 0;
EvtVectorParticle parent{};
parent.noLifeTime();
parent.init( getParentId(),
EvtVector4R( EvtPDL::getMass( getParentId() ), 0, 0, 0 ) );
parent.setDiagonalSpinDensity();
EvtAmp amp;
EvtId daughters[3] = { getDaug( 0 ), getDaug( 1 ), getDaug( 2 ) };
amp.init( getParentId(), 3, daughters );
parent.makeDaughters( 3, daughters );
EvtParticle* scalar = parent.getDaug( 0 );
EvtParticle* lep1 = parent.getDaug( 1 );
EvtParticle* lep2 = parent.getDaug( 2 );
scalar->noLifeTime();
lep1->noLifeTime();
lep2->noLifeTime();
EvtSpinDensity rho;
rho.setDiag( parent.getSpinStates() );
const double M0 = EvtPDL::getMass( getParentId() );
const double mL = EvtPDL::getMass( getDaug( 0 ) );
const double m1 = EvtPDL::getMass( getDaug( 1 ) );
const double m2 = EvtPDL::getMass( getDaug( 2 ) );
const double m12Sum = m1 + m2;
const double m12Del = m1 - m2;
const double q2min = ( m1 + m2 ) * ( m1 + m2 );
const double q2max = ( M0 - mL ) * ( M0 - mL );
const double mSqSum = M0 * M0 + mL * mL;
EvtVector4R p4scalar, p4lep1, p4lep2, boost;
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< " EvtVtoSll is probing whole phase space ..." << std::endl;
double prob = 0;
const int nsteps = 5000;
for ( int i = 0; i <= nsteps; i++ ) {
const double q2 = q2min + i * ( q2max - q2min ) / nsteps;
const double eScalar = ( mSqSum - q2 ) / ( 2 * M0 );
const double pstar = i == 0 ? 0
: sqrt( q2 - m12Sum * m12Sum ) *
sqrt( q2 - m12Del * m12Del ) /
( 2 * sqrt( q2 ) );
boost.set( M0 - eScalar, 0, 0, +sqrt( eScalar * eScalar - mL * mL ) );
if ( i != nsteps ) {
p4scalar.set( eScalar, 0, 0, -sqrt( eScalar * eScalar - mL * mL ) );
} else {
p4scalar.set( mL, 0, 0, 0 );
} }
void EvtVtoSll::decay( EvtParticle* p ) const double pstarSq = pstar * pstar;
const double E1star = sqrt( pstarSq + m1 * m1 );
const double E2star = sqrt( pstarSq + m2 * m2 );
for ( int j = 0; j <= 45; j++ ) {
const double theta = j * EvtConst::pi / 45;
const double pstarT = pstar * sin( theta );
const double pstarZ = pstar * cos( theta );
p4lep1.set( E1star, 0, pstarT, pstarZ );
p4lep2.set( E2star, 0, -pstarT, -pstarZ );
if ( i != nsteps ) // At maximal q2 we are already in correct frame as scalar child and W/Zvirtual are at rest
{ {
p->initializePhaseSpace( getNDaug(), getDaugs() ); p4lep1 = boostTo( p4lep1, boost );
p4lep2 = boostTo( p4lep2, boost );
}
scalar->init( getDaug( 0 ), p4scalar );
lep1->init( getDaug( 1 ), p4lep1 );
lep2->init( getDaug( 2 ), p4lep2 );
calcAmp( parent, amp );
prob = rho.normalizedProb( amp.getSpinDensity() );
// In case of electron mode add pole
if ( m_electronMode ) {
prob /= 1.0 + m_poleSize / ( q2 * q2 );
}
if ( prob > theProbMax ) {
theProbMax = prob;
theProbMax_q2 = q2;
theProbMax_theta = theta;
tlathamUnsubmitted
Not Done
Inline Actions

I've just run this revision through the CI and I notice that this produces quite a lot of printout. Can we put this inside some sort of if ( m_debug ) condition so that it isn't normally printed?

tlatham: I've just run this revision through the CI and I notice that this produces quite a lot of…
}
}
}
theProbMax *= 1.01;
setProbMax( theProbMax );
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< " EvtVtoSll set up maximum probability to " << theProbMax
<< " found at q2 = " << theProbMax_q2 << " ("
<< nsteps * ( theProbMax_q2 - q2min ) / ( q2max - q2min )
<< " %) and theta = " << theProbMax_theta * 180 / EvtConst::pi
<< std::endl;
}
EvtParticle *l1, *l2; void EvtVtoSll::decay( EvtParticle* parent )
l1 = p->getDaug( 1 ); {
l2 = p->getDaug( 2 ); // Phase space initialization depends on what leptons are
if ( m_electronMode ) {
setWeight( parent->initializePhaseSpace( getNDaug(), getDaugs(), false,
m_poleSize, 1, 2 ) );
} else {
parent->initializePhaseSpace( getNDaug(), getDaugs() );
}
EvtVector4C l11, l12, l21, l22; calcAmp( *parent, _amp2 );
l11 = EvtLeptonVCurrent( l1->spParent( 0 ), l2->spParent( 0 ) ); }
l12 = EvtLeptonVCurrent( l1->spParent( 0 ), l2->spParent( 1 ) );
l21 = EvtLeptonVCurrent( l1->spParent( 1 ), l2->spParent( 0 ) );
l22 = EvtLeptonVCurrent( l1->spParent( 1 ), l2->spParent( 1 ) );
EvtVector4C eps0 = p->eps( 0 ); void EvtVtoSll::calcAmp( const EvtParticle& parent, EvtAmp& amp ) const
EvtVector4C eps1 = p->eps( 1 ); {
EvtVector4C eps2 = p->eps( 2 ); const EvtParticle& l1 = *parent.getDaug( 1 );
const EvtParticle& l2 = *parent.getDaug( 2 );
EvtVector4R P = p->getP4Restframe(); const EvtVector4C l11 = EvtLeptonVCurrent( l1.spParent( 0 ),
EvtVector4R k = l1->getP4() + l2->getP4(); l2.spParent( 0 ) );
double k2 = k * k; const EvtVector4C l12 = EvtLeptonVCurrent( l1.spParent( 0 ),
l2.spParent( 1 ) );
const EvtVector4C l21 = EvtLeptonVCurrent( l1.spParent( 1 ),
l2.spParent( 0 ) );
const EvtVector4C l22 = EvtLeptonVCurrent( l1.spParent( 1 ),
l2.spParent( 1 ) );
const EvtVector4C eps0 = parent.eps( 0 );
const EvtVector4C eps1 = parent.eps( 1 );
const EvtVector4C eps2 = parent.eps( 2 );
const EvtVector4R P = parent.getP4Restframe();
const EvtVector4R k = l1.getP4() + l2.getP4();
const double k2 = k * k;
EvtTensor4C T( dual( EvtGenFunctions::directProd( P, ( 1.0 / k2 ) * k ) ) ); const EvtTensor4C T(
dual( EvtGenFunctions::directProd( P, ( 1.0 / k2 ) * k ) ) );
double M2 = p->mass(); double M2 = parent.mass();
M2 *= M2; M2 *= M2;
double m2 = l1->mass(); double m2 = l1.mass();
m2 *= m2; m2 *= m2;
double norm = 1.0 / sqrt( 2 * M2 + 4 * m2 - 4 * m2 * m2 / M2 ); const double norm = 1.0 / sqrt( 2 * M2 + 4 * m2 - 4 * m2 * m2 / M2 );
vertex( 0, 0, 0, norm * ( eps0 * T.cont2( l11 ) ) ); amp.vertex( 0, 0, 0, norm * ( eps0 * T.cont2( l11 ) ) );
vertex( 0, 0, 1, norm * ( eps0 * T.cont2( l12 ) ) ); amp.vertex( 0, 0, 1, norm * ( eps0 * T.cont2( l12 ) ) );
vertex( 0, 1, 0, norm * ( eps0 * T.cont2( l21 ) ) ); amp.vertex( 0, 1, 0, norm * ( eps0 * T.cont2( l21 ) ) );
vertex( 0, 1, 1, norm * ( eps0 * T.cont2( l22 ) ) ); amp.vertex( 0, 1, 1, norm * ( eps0 * T.cont2( l22 ) ) );
vertex( 1, 0, 0, norm * ( eps1 * T.cont2( l11 ) ) ); amp.vertex( 1, 0, 0, norm * ( eps1 * T.cont2( l11 ) ) );
vertex( 1, 0, 1, norm * ( eps1 * T.cont2( l12 ) ) ); amp.vertex( 1, 0, 1, norm * ( eps1 * T.cont2( l12 ) ) );
vertex( 1, 1, 0, norm * ( eps1 * T.cont2( l21 ) ) ); amp.vertex( 1, 1, 0, norm * ( eps1 * T.cont2( l21 ) ) );
vertex( 1, 1, 1, norm * ( eps1 * T.cont2( l22 ) ) ); amp.vertex( 1, 1, 1, norm * ( eps1 * T.cont2( l22 ) ) );
vertex( 2, 0, 0, norm * ( eps2 * T.cont2( l11 ) ) ); amp.vertex( 2, 0, 0, norm * ( eps2 * T.cont2( l11 ) ) );
vertex( 2, 0, 1, norm * ( eps2 * T.cont2( l12 ) ) ); amp.vertex( 2, 0, 1, norm * ( eps2 * T.cont2( l12 ) ) );
vertex( 2, 1, 0, norm * ( eps2 * T.cont2( l21 ) ) ); amp.vertex( 2, 1, 0, norm * ( eps2 * T.cont2( l21 ) ) );
vertex( 2, 1, 1, norm * ( eps2 * T.cont2( l22 ) ) ); amp.vertex( 2, 1, 1, norm * ( eps2 * T.cont2( l22 ) ) );
return; return;
} }