Page MenuHomeHEPForge
Differential D58 Diff 259 src/EvtGenBase/EvtAbsLineShape.cpp

Changeset View

Standalone View

View Options

src/EvtGenBase/EvtAbsLineShape.cpp

Show All 32 Lines
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
using namespace std; using namespace std;
EvtAbsLineShape::EvtAbsLineShape( double mass, double width, double maxRange, EvtAbsLineShape::EvtAbsLineShape( double mass, double width, double maxRange,
EvtSpinType::spintype sp ) EvtSpinType::spintype sp )
{ {
_includeDecayFact = false; m_includeDecayFact = false;
_includeBirthFact = false; m_includeBirthFact = false;
_mass = mass; m_mass = mass;
_width = width; m_width = width;
_spin = sp; m_spin = sp;
_maxRange = maxRange; m_maxRange = maxRange;
double maxdelta = 15.0 * width; double maxdelta = 15.0 * width;
//if ( width>0.001 ) { //if ( width>0.001 ) {
// if ( 5.0*width < 0.6 ) maxdelta = 0.6; // if ( 5.0*width < 0.6 ) maxdelta = 0.6;
//} //}
if ( maxRange > 0.00001 ) { if ( maxRange > 0.00001 ) {
_massMax = mass + maxdelta; m_massMax = mass + maxdelta;
_massMin = mass - maxRange; m_massMin = mass - maxRange;
} else { } else {
_massMax = mass + maxdelta; m_massMax = mass + maxdelta;
_massMin = mass - 15.0 * width; m_massMin = mass - 15.0 * width;
} }
if ( _massMin < 0. ) if ( m_massMin < 0. )
_massMin = 0.; m_massMin = 0.;
_massMax = mass + maxdelta; m_massMax = mass + maxdelta;
} }
EvtAbsLineShape::EvtAbsLineShape( const EvtAbsLineShape& x ) EvtAbsLineShape::EvtAbsLineShape( const EvtAbsLineShape& x )
{ {
_includeDecayFact = x._includeDecayFact; m_includeDecayFact = x.m_includeDecayFact;
_includeBirthFact = x._includeBirthFact; m_includeBirthFact = x.m_includeBirthFact;
_mass = x._mass; m_mass = x.m_mass;
_massMax = x._massMax; m_massMax = x.m_massMax;
_massMin = x._massMin; m_massMin = x.m_massMin;
_width = x._width; m_width = x.m_width;
_spin = x._spin; m_spin = x.m_spin;
_maxRange = x._maxRange; m_maxRange = x.m_maxRange;
} }
EvtAbsLineShape& EvtAbsLineShape::operator=( const EvtAbsLineShape& x ) EvtAbsLineShape& EvtAbsLineShape::operator=( const EvtAbsLineShape& x )
{ {
_includeDecayFact = x._includeDecayFact; m_includeDecayFact = x.m_includeDecayFact;
_includeBirthFact = x._includeBirthFact; m_includeBirthFact = x.m_includeBirthFact;
_mass = x._mass; m_mass = x.m_mass;
_massMax = x._massMax; m_massMax = x.m_massMax;
_massMin = x._massMin; m_massMin = x.m_massMin;
_width = x._width; m_width = x.m_width;
_spin = x._spin; m_spin = x.m_spin;
_maxRange = x._maxRange; m_maxRange = x.m_maxRange;
return *this; return *this;
} }
EvtAbsLineShape* EvtAbsLineShape::clone() EvtAbsLineShape* EvtAbsLineShape::clone()
{ {
return new EvtAbsLineShape( *this ); return new EvtAbsLineShape( *this );
} }
double EvtAbsLineShape::rollMass() double EvtAbsLineShape::rollMass()
{ {
double ymin, ymax; double ymin, ymax;
double temp; double temp;
if ( _width < 0.0001 ) { if ( m_width < 0.0001 ) {
return _mass; return m_mass;
} else { } else {
ymin = atan( 2.0 * ( _massMin - _mass ) / _width ); ymin = atan( 2.0 * ( m_massMin - m_mass ) / m_width );
ymax = atan( 2.0 * ( _massMax - _mass ) / _width ); ymax = atan( 2.0 * ( m_massMax - m_mass ) / m_width );
temp = ( _mass + temp = ( m_mass +
( ( _width / 2.0 ) * tan( EvtRandom::Flat( ymin, ymax ) ) ) ); ( ( m_width / 2.0 ) * tan( EvtRandom::Flat( ymin, ymax ) ) ) );
return temp; return temp;
} }
} }
double EvtAbsLineShape::getRandMass( EvtId* parId, int /* nDaug */, double EvtAbsLineShape::getRandMass( EvtId* parId, int /* nDaug */,
EvtId* /*dauId*/, EvtId* /*othDaugId*/, EvtId* /*dauId*/, EvtId* /*othDaugId*/,
double maxMass, double* /*dauMasses*/ ) double maxMass, double* /*dauMasses*/ )
{ {
if ( _width < 0.0001 ) if ( m_width < 0.0001 )
return _mass; return m_mass;
//its not flat - but generated according to a BW //its not flat - but generated according to a BW
if ( maxMass > 0 && maxMass < _massMin ) { if ( maxMass > 0 && maxMass < m_massMin ) {
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
<< "In EvtAbsLineShape::getRandMass:" << endl; << "In EvtAbsLineShape::getRandMass:" << endl;
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
<< "Cannot create a particle with a minimal mass of " << _massMin << "Cannot create a particle with a minimal mass of " << m_massMin
<< " from a " << EvtPDL::name( *parId ) << " from a " << EvtPDL::name( *parId )
<< " decay with available left-over mass-energy " << maxMass << " decay with available left-over mass-energy " << maxMass
<< ". Returning 0.0 mass. The rest of this decay chain will probably fail..." << ". Returning 0.0 mass. The rest of this decay chain will probably fail..."
<< endl; << endl;
return 0.0; return 0.0;
} }
double mMin = _massMin; double mMin = m_massMin;
double mMax = _massMax; double mMax = m_massMax;
if ( maxMass > -0.5 && maxMass < mMax ) if ( maxMass > -0.5 && maxMass < mMax )
mMax = maxMass; mMax = maxMass;
double ymin = atan( 2.0 * ( mMin - _mass ) / _width ); double ymin = atan( 2.0 * ( mMin - m_mass ) / m_width );
double ymax = atan( 2.0 * ( mMax - _mass ) / _width ); double ymax = atan( 2.0 * ( mMax - m_mass ) / m_width );
return ( _mass + ( ( _width / 2.0 ) * tan( EvtRandom::Flat( ymin, ymax ) ) ) ); return ( m_mass +
// return EvtRandom::Flat(_massMin,_massMax); ( ( m_width / 2.0 ) * tan( EvtRandom::Flat( ymin, ymax ) ) ) );
// return EvtRandom::Flat(m_massMin,m_massMax);
} }
double EvtAbsLineShape::getMassProb( double mass, double massPar, int nDaug, double EvtAbsLineShape::getMassProb( double mass, double massPar, int nDaug,
double* massDau ) double* massDau )
{ {
double dTotMass = 0.; double dTotMass = 0.;
if ( nDaug > 1 ) { if ( nDaug > 1 ) {
int i; int i;
for ( i = 0; i < nDaug; i++ ) { for ( i = 0; i < nDaug; i++ ) {
dTotMass += massDau[i]; dTotMass += massDau[i];
} }
//EvtGenReport(EVTGEN_INFO,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl; //EvtGenReport(EVTGEN_INFO,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl;
// if ( (mass-dTotMass)<0.0001 ) return 0.; // if ( (mass-dTotMass)<0.0001 ) return 0.;
if ( ( mass < dTotMass ) ) if ( ( mass < dTotMass ) )
return 0.; return 0.;
} }
if ( _width < 0.0001 ) if ( m_width < 0.0001 )
return 1.; return 1.;
// no parent - lets not panic // no parent - lets not panic
if ( massPar > 0.0000000001 ) { if ( massPar > 0.0000000001 ) {
if ( mass > massPar ) if ( mass > massPar )
return 0.; return 0.;
} }
//Otherwise return the right value. //Otherwise return the right value.
//Fortunately we have generated events according to a non-rel BW, so //Fortunately we have generated events according to a non-rel BW, so
//just return.. //just return..
//EvtPropBreitWigner bw(_mass,_width); //EvtPropBreitWigner bw(m_mass,m_width);
//EvtPropFactor<EvtTwoBodyVertex> f(bw); //EvtPropFactor<EvtTwoBodyVertex> f(bw);
//EvtComplex fm=f.eval(mass); //EvtComplex fm=f.eval(mass);
//EvtComplex fm0=f.eval(_mass); //EvtComplex fm0=f.eval(m_mass);
//return (abs(fm)*abs(fm))/(abs(fm0)*abs(fm0)); //return (abs(fm)*abs(fm))/(abs(fm0)*abs(fm0));
return 1.0; return 1.0;
} }