Differential D98 Diff 424 src/EvtGenModels/EvtPropSLPole.cpp

Changeset View

Standalone View

src/EvtGenModels/EvtPropSLPole.cpp

Show First 20 Lines • Show All 197 Lines • ▼ Show 20 Lines	double EvtPropSLPole::calBreitWigner( EvtParticle* pmeson, EvtPoint1D point )
if ( _massMin < 0. )		if ( _massMin < 0. )
_massMin = 0.;		_massMin = 0.;

EvtParticle* par = pmeson->getParent();		EvtParticle* par = pmeson->getParent();
double maxMass = -1.;		double maxMass = -1.;
if ( par != nullptr ) {		if ( par != nullptr ) {
if ( par->hasValidP4() )		if ( par->hasValidP4() )
maxMass = par->mass();		maxMass = par->mass();
for ( size_t i = 0; i < par->getNDaug(); i++ ) {		for ( std::size_t i = 0; i < par->getNDaug(); i++ ) {
EvtParticle* tDaug = par->getDaug( i );		EvtParticle* tDaug = par->getDaug( i );
if ( pmeson != tDaug )		if ( pmeson != tDaug )
maxMass -= EvtPDL::getMinMass( tDaug->getId() );		maxMass -= EvtPDL::getMinMass( tDaug->getId() );
}		}
}		}

EvtId* dauId = nullptr;		std::vector<EvtId> dauId;
double* dauMasses = nullptr;		std::vector<double> dauMasses;
size_t nDaug = pmeson->getNDaug();		const std::size_t nDaug{ pmeson->getNDaug() };
if ( nDaug > 0 ) {		if ( nDaug > 0 ) {
dauId = new EvtId[nDaug];		dauId.resize( nDaug );
dauMasses = new double[nDaug];		dauMasses.resize( nDaug );
for ( size_t j = 0; j < nDaug; j++ ) {		for ( std::size_t j = 0; j < nDaug; j++ ) {
dauId[j] = pmeson->getDaug( j )->getId();		dauId[j] = pmeson->getDaug( j )->getId();
dauMasses[j] = pmeson->getDaug( j )->mass();		dauMasses[j] = pmeson->getDaug( j )->mass();
}		}
}		}
EvtId* parId = nullptr;		EvtId parId;
EvtId* othDaugId = nullptr;		EvtId othDaugId;
EvtParticle* tempPar = pmeson->getParent();		EvtParticle* tempPar = pmeson->getParent();
if ( tempPar ) {		if ( tempPar ) {
parId = new EvtId( tempPar->getId() );		parId = tempPar->getId();
if ( tempPar->getNDaug() == 2 ) {		if ( tempPar->getNDaug() == 2 ) {
if ( tempPar->getDaug( 0 ) == pmeson )		if ( tempPar->getDaug( 0 ) == pmeson ) {
othDaugId = new EvtId( tempPar->getDaug( 1 )->getId() );		othDaugId = tempPar->getDaug( 1 )->getId();
else		} else {
othDaugId = new EvtId( tempPar->getDaug( 0 )->getId() );		othDaugId = tempPar->getDaug( 0 )->getId();
		}
}		}
}		}

if ( nDaug != 2 )		if ( nDaug != 2 )
return calBreitWignerBasic( maxMass );		return calBreitWignerBasic( maxMass );

if ( _width < 0.00001 )		if ( _width < 0.00001 )
return 1.0;		return 1.0;
Show All 34 Lines	double EvtPropSLPole::calBreitWigner( EvtParticle* pmeson, EvtPoint1D point )
// I'm not sure how to define the vertex factor here - so retreat to nonRel code.		// I'm not sure how to define the vertex factor here - so retreat to nonRel code.
if ( ( massD1 + massD2 ) > _mass )		if ( ( massD1 + massD2 ) > _mass )
return calBreitWignerBasic( maxMass );		return calBreitWignerBasic( maxMass );

//parent vertex factor not yet implemented		//parent vertex factor not yet implemented
double massOthD = -10.;		double massOthD = -10.;
double massParent = -10.;		double massParent = -10.;
int birthl = -10;		int birthl = -10;
if ( othDaugId ) {		if ( othDaugId != EvtId{} ) {
EvtSpinType::spintype spinOth = EvtPDL::getSpinType( *othDaugId );		EvtSpinType::spintype spinOth = EvtPDL::getSpinType( othDaugId );
EvtSpinType::spintype spinPar = EvtPDL::getSpinType( *parId );		EvtSpinType::spintype spinPar = EvtPDL::getSpinType( parId );

int tt1 = EvtSpinType::getSpin2( spinOth );		int tt1 = EvtSpinType::getSpin2( spinOth );
int tt2 = EvtSpinType::getSpin2( spinPar );		int tt2 = EvtSpinType::getSpin2( spinPar );
int tt3 = EvtSpinType::getSpin2( _spin );		int tt3 = EvtSpinType::getSpin2( _spin );

//figure the min and max allowwed "spins" for the daughters state		//figure the min and max allowwed "spins" for the daughters state
if ( ( tt1 <= 4 ) && ( tt2 <= 4 ) ) {		if ( ( tt1 <= 4 ) && ( tt2 <= 4 ) ) {
birthl = std::max( tt3 - tt2 - tt1,		birthl = std::max( tt3 - tt2 - tt1,
std::max( tt2 - tt3 - tt1, tt1 - tt3 - tt2 ) );		std::max( tt2 - tt3 - tt1, tt1 - tt3 - tt2 ) );
if ( birthl < 0 )		if ( birthl < 0 )
birthl = 0;		birthl = 0;

massOthD = EvtPDL::getMeanMass( *othDaugId );		massOthD = EvtPDL::getMeanMass( othDaugId );
massParent = EvtPDL::getMeanMass( *parId );		massParent = EvtPDL::getMeanMass( parId );
}		}
}		}
double massM = _massMax;		double massM = _massMax;
if ( ( maxMass > -0.5 ) && ( maxMass < massM ) )		if ( ( maxMass > -0.5 ) && ( maxMass < massM ) )
massM = maxMass;		massM = maxMass;

//special case... if the parent mass is _fixed_ we can do a little better		//special case... if the parent mass is _fixed_ we can do a little better
//and only for a two body decay as that seems to be where we have problems		//and only for a two body decay as that seems to be where we have problems
Show All 18 Lines	if ( massParent > -1. ) {
amp.addBirthFactFF();		amp.addBirthFactFF();
}		}
}		}

EvtAmpPdf<EvtPoint1D> pdf( amp );		EvtAmpPdf<EvtPoint1D> pdf( amp );

double ampVal = sqrt( pdf.evaluate( point ) );		double ampVal = sqrt( pdf.evaluate( point ) );

if ( parId )
delete parId;
if ( othDaugId )
delete othDaugId;
if ( dauId )
delete[] dauId;
if ( dauMasses )
delete[] dauMasses;

return ampVal;		return ampVal;
}		}

double EvtPropSLPole::calcMaxProb( EvtId parent, EvtId meson, EvtId lepton,		double EvtPropSLPole::calcMaxProb( EvtId parent, EvtId meson, EvtId lepton,
EvtId nudaug, EvtSemiLeptonicFF* FormFactors )		EvtId nudaug, EvtSemiLeptonicFF* FormFactors )
{		{
//This routine takes the arguements parent, meson, and lepton		//This routine takes the arguements parent, meson, and lepton
//number, and a form factor model, and returns a maximum		//number, and a form factor model, and returns a maximum
Show All 13 Lines	double EvtPropSLPole::calcMaxProb( EvtId parent, EvtId meson, EvtId lepton,

//cludge to avoid generating random numbers!		//cludge to avoid generating random numbers!
scalar_part->noLifeTime();		scalar_part->noLifeTime();

EvtVector4R p_init;		EvtVector4R p_init;

p_init.set( EvtPDL::getMass( parent ), 0.0, 0.0, 0.0 );		p_init.set( EvtPDL::getMass( parent ), 0.0, 0.0, 0.0 );
scalar_part->init( parent, p_init );		scalar_part->init( parent, p_init );
root_part = (EvtParticle*)scalar_part;		root_part = scalar_part;
// root_part->set_type(EvtSpinType::SCALAR);		// root_part->set_type(EvtSpinType::SCALAR);
root_part->setDiagonalSpinDensity();		root_part->setDiagonalSpinDensity();

EvtParticle daughter, lep, *trino;		EvtParticle daughter, lep, *trino;

EvtAmp amp;		EvtAmp amp;

EvtId listdaug[3];		EvtId listdaug[3];
▲ Show 20 Lines • Show All 94 Lines • ▼ Show 20 Lines	for ( massiter = 0; massiter < 3; massiter++ ) {
//Now initialize the daughters...		//Now initialize the daughters...

daughter->init( meson, p4meson );		daughter->init( meson, p4meson );
lep->init( lepton, p4lepton );		lep->init( lepton, p4lepton );
trino->init( nudaug, p4nu );		trino->init( nudaug, p4nu );

calcamp->CalcAmp( root_part, amp, FormFactors );		calcamp->CalcAmp( root_part, amp, FormFactors );

EvtPoint1D* point = new EvtPoint1D( mass[0] );		EvtPoint1D point( mass[0] );

double meson_BWAmp = calBreitWigner( daughter, *point );		double meson_BWAmp = calBreitWigner( daughter, point );

int list[2];		int list[2];
list[0] = 0;		list[0] = 0;
list[1] = 0;		list[1] = 0;
amp.vertex( 0, 0, amp.getAmp( list ) * meson_BWAmp );		amp.vertex( 0, 0, amp.getAmp( list ) * meson_BWAmp );
list[0] = 0;		list[0] = 0;
list[1] = 1;		list[1] = 1;
amp.vertex( 0, 1, amp.getAmp( list ) * meson_BWAmp );		amp.vertex( 0, 1, amp.getAmp( list ) * meson_BWAmp );
▲ Show 20 Lines • Show All 65 Lines • Show Last 20 Lines