Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenBase/EvtAmp.cpp
Show First 20 Lines • Show All 135 Lines • ▼ Show 20 Lines | for ( int i = 1; i < _nontrivial; i++ ) { | ||||
position += nstatepad * ind[i]; | position += nstatepad * ind[i]; | ||||
} | } | ||||
return _amp[position]; | return _amp[position]; | ||||
} | } | ||||
EvtSpinDensity EvtAmp::getSpinDensity() | EvtSpinDensity EvtAmp::getSpinDensity() | ||||
{ | { | ||||
int np = _pstates; | |||||
EvtSpinDensity rho; | EvtSpinDensity rho; | ||||
rho.setDim( _pstates ); | rho.setDim( np ); | ||||
if ( np == 1 ) { | |||||
EvtComplex temp; | double a = abs2( _amp[0] ); | ||||
int n = 1; | |||||
int i, j, n; | for ( int i = 0; i < _nontrivial; i++ ) | ||||
if ( _pstates == 1 ) { | |||||
if ( _nontrivial == 0 ) { | |||||
rho.set( 0, 0, _amp[0] * conj( _amp[0] ) ); | |||||
return rho; | |||||
} | |||||
n = 1; | |||||
temp = EvtComplex( 0.0 ); | |||||
for ( i = 0; i < _nontrivial; i++ ) { | |||||
n *= _nstate[i]; | n *= _nstate[i]; | ||||
for ( int i = 1; i < n; i++ ) | |||||
a += abs2( _amp[i] ); | |||||
rho.set( 0, 0, EvtComplex( a ) ); | |||||
} else { | |||||
int n = 1; | |||||
for ( int i = 0; i < _ndaug; i++ ) | |||||
n *= dstates[i]; | |||||
for ( int i = 0; i < np; i++ ) { | |||||
for ( int j = 0; j < np; j++ ) { | |||||
EvtComplex a; | |||||
for ( int k = 0; k < n; k++ ) | |||||
a += _amp[np * k + i] * conj( _amp[np * k + j] ); | |||||
rho.set( i, j, a ); | |||||
} | } | ||||
for ( i = 0; i < n; i++ ) { | |||||
temp += _amp[i] * conj( _amp[i] ); | |||||
} | |||||
rho.set( 0, 0, temp ); | |||||
; | |||||
return rho; | |||||
} | |||||
else { | |||||
for ( i = 0; i < _pstates; i++ ) { | |||||
for ( j = 0; j < _pstates; j++ ) { | |||||
temp = EvtComplex( 0.0 ); | |||||
int kk; | |||||
int allloop = 1; | |||||
for ( kk = 0; kk < _ndaug; kk++ ) { | |||||
allloop *= dstates[kk]; | |||||
} | |||||
for ( kk = 0; kk < allloop; kk++ ) { | |||||
temp += _amp[_pstates * kk + i] * | |||||
conj( _amp[_pstates * kk + j] ); | |||||
} | |||||
// if (_nontrivial>3){ | |||||
//EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Can't handle so many states in EvtAmp!"<<endl; | |||||
//} | |||||
rho.set( i, j, temp ); | |||||
} | } | ||||
} | } | ||||
return rho; | return rho; | ||||
} | } | ||||
} | |||||
EvtSpinDensity EvtAmp::getBackwardSpinDensity( EvtSpinDensity* rho_list ) | EvtSpinDensity EvtAmp::getBackwardSpinDensity( EvtSpinDensity* rho_list ) | ||||
{ | { | ||||
EvtSpinDensity rho; | EvtSpinDensity rho; | ||||
rho.setDim( _pstates ); | rho.setDim( _pstates ); | ||||
if ( _pstates == 1 ) { | if ( _pstates == 1 ) { | ||||
▲ Show 20 Lines • Show All 42 Lines • ▼ Show 20 Lines | for ( k = 0; k < i; k++ ) { | ||||
if ( dstates[k] != 1 ) { | if ( dstates[k] != 1 ) { | ||||
ampprime = ampprime.contract( _dnontrivial[k], rho_list[k + 1] ); | ampprime = ampprime.contract( _dnontrivial[k], rho_list[k + 1] ); | ||||
} | } | ||||
} | } | ||||
return ampprime.contract( _dnontrivial[i], ( *this ) ); | return ampprime.contract( _dnontrivial[i], ( *this ) ); | ||||
} | } | ||||
EvtAmp EvtAmp::contract( int k, const EvtSpinDensity& rho ) | EvtAmp EvtAmp::contract( int istate, const EvtSpinDensity& rho ) | ||||
{ | { | ||||
// contract amplitude tensor with spin density tensor | |||||
// istate -- contraction index | |||||
// dimension of rho has to be _nstate[istate] | |||||
EvtAmp temp; | EvtAmp temp; | ||||
int i, j; | |||||
temp._ndaug = _ndaug; | temp._ndaug = _ndaug; | ||||
temp._pstates = _pstates; | temp._pstates = _pstates; | ||||
temp._nontrivial = _nontrivial; | temp._nontrivial = _nontrivial; | ||||
for ( int i = 0; i < _ndaug; i++ ) { | |||||
for ( i = 0; i < _ndaug; i++ ) { | |||||
temp.dstates[i] = dstates[i]; | temp.dstates[i] = dstates[i]; | ||||
temp._dnontrivial[i] = _dnontrivial[i]; | temp._dnontrivial[i] = _dnontrivial[i]; | ||||
} | } | ||||
for ( int i = 0; i < _nontrivial; i++ ) | |||||
if ( _nontrivial == 0 ) { | |||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | |||||
<< "Should not be here EvtAmp!" << endl; | |||||
} | |||||
for ( i = 0; i < _nontrivial; i++ ) { | |||||
temp._nstate[i] = _nstate[i]; | temp._nstate[i] = _nstate[i]; | ||||
} | int nup = 1, stride = 1; | ||||
for ( int i = istate + 1; i < _nontrivial; i++ ) | |||||
nup *= _nstate[i]; | |||||
for ( int i = 0; i < istate; i++ ) | |||||
stride *= _nstate[i]; | |||||
int ns = _nstate[istate], offstride = stride * ns, namp = offstride * nup; | |||||
#if 1 | |||||
for ( int i = 0; i < stride; i++ ) { | |||||
for ( int off = i; off < namp; off += offstride ) { | |||||
for ( int k = 0; k < ns; k++ ) { | |||||
EvtComplex c; | EvtComplex c; | ||||
for ( int j = 0; j < ns; j++ ) | |||||
int index[10]; | c += rho( j, k ) * _amp[off + j * stride]; | ||||
for ( i = 0; i < 10; i++ ) { | temp._amp[off + k * stride] = c; | ||||
index[i] = 0; | |||||
} | } | ||||
int allloop = 1; | |||||
int indflag, ii; | |||||
for ( i = 0; i < _nontrivial; i++ ) { | |||||
allloop *= _nstate[i]; | |||||
} | } | ||||
} | |||||
for ( i = 0; i < allloop; i++ ) { | // temp.dump(); | ||||
c = EvtComplex( 0.0 ); | #endif | ||||
int tempint = index[k]; | #if 0 | ||||
for ( j = 0; j < _nstate[k]; j++ ) { | int index[10] = {0}; | ||||
index[k] = j; | for (int i = 0; i < namp; i++ ) { | ||||
EvtComplex c; | |||||
int tempint = index[istate]; | |||||
for (int j = 0; j < _nstate[istate]; j++ ) { | |||||
// std::cout<<i<<" "<<j<<" "<<tempint<<std::endl; | |||||
index[istate] = j; | |||||
c += rho.get( j, tempint ) * getAmp( index ); | c += rho.get( j, tempint ) * getAmp( index ); | ||||
} | } | ||||
index[k] = tempint; | index[istate] = tempint; | ||||
temp.setAmp( index, c ); | temp.setAmp( index, c ); | ||||
int indflag = 0; | |||||
indflag = 0; | for (int j = 0; j < _nontrivial; j++ ) { | ||||
for ( ii = 0; ii < _nontrivial; ii++ ) { | |||||
if ( indflag == 0 ) { | if ( indflag == 0 ) { | ||||
if ( index[ii] == ( _nstate[ii] - 1 ) ) { | if ( index[j] == ( _nstate[j] - 1 ) ) { | ||||
index[ii] = 0; | index[j] = 0; | ||||
} else { | } else { | ||||
indflag = 1; | indflag = 1; | ||||
index[ii] += 1; | index[j] += 1; | ||||
} | } | ||||
} | } | ||||
} | } | ||||
} | } | ||||
// temp.dump(); | |||||
#endif | |||||
return temp; | return temp; | ||||
} | } | ||||
EvtSpinDensity EvtAmp::contract( int k, const EvtAmp& amp2 ) | EvtSpinDensity EvtAmp::contract( int istate, const EvtAmp& amp2 ) | ||||
{ | { | ||||
#if 1 | |||||
int nup = 1, stride = 1; | |||||
for ( int i = istate + 1; i < _nontrivial; i++ ) | |||||
nup *= _nstate[i]; | |||||
for ( int i = 0; i < istate; i++ ) | |||||
stride *= _nstate[i]; | |||||
int ns = _nstate[istate], offstride = stride * ns, namp = offstride * nup; | |||||
EvtSpinDensity mrho( ns ); | |||||
for ( int i = 0; i < stride; i++ ) { | |||||
for ( int off = i; off < namp; off += offstride ) { | |||||
for ( int k = 0; k < ns; k++ ) { | |||||
EvtComplex ak = conj( amp2._amp[off + k * stride] ); | |||||
for ( int j = k; j < ns; j++ ) { | |||||
EvtComplex aj = _amp[off + j * stride]; | |||||
mrho( j, k ) += aj * ak; | |||||
// if(k==j)std::cout<<" "<<i<<" "<<off<<" "<<k<<" "<<ak<<" "<<aj<<" "<<aj*ak<<std::endl; | |||||
} | |||||
} | |||||
} | |||||
} | |||||
for ( int k = 0; k < ns; k++ ) { | |||||
for ( int j = 0; j < k; j++ ) { | |||||
mrho( j, k ) = conj( mrho( k, j ) ); | |||||
} | |||||
} | |||||
// std::cout<<mrho; | |||||
return mrho; | |||||
#endif | |||||
#if 0 | |||||
int k = istate; | |||||
int i, j, l; | int i, j, l; | ||||
EvtComplex temp; | EvtComplex temp; | ||||
EvtSpinDensity rho; | EvtSpinDensity rho( _nstate[k] ); | ||||
// rho.setDim( _nstate[k] ); | |||||
rho.setDim( _nstate[k] ); | |||||
int allloop = 1; | int allloop = 1; | ||||
int indflag, ii; | int indflag, ii; | ||||
for ( i = 0; i < _nontrivial; i++ ) { | for ( i = 0; i < _nontrivial; i++ ) { | ||||
allloop *= _nstate[i]; | allloop *= _nstate[i]; | ||||
} | } | ||||
int index[10]; | int index[10]; | ||||
Show All 33 Lines | for ( i = 0; i < _nstate[k]; i++ ) { | ||||
} | } | ||||
} | } | ||||
} | } | ||||
} | } | ||||
} | } | ||||
rho.set( i, j, temp ); | rho.set( i, j, temp ); | ||||
} | } | ||||
} | } | ||||
// std::cout<<rho-mrho; | |||||
return rho; | return rho; | ||||
#endif | |||||
} | } | ||||
EvtAmp EvtAmp::contract( int, const EvtAmp&, const EvtAmp& ) | EvtAmp EvtAmp::contract( int, const EvtAmp&, const EvtAmp& ) | ||||
{ | { | ||||
//Do we need this method? | //Do we need this method? | ||||
EvtAmp tmp; | EvtAmp tmp; | ||||
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) | EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) | ||||
<< "EvtAmp::contract not written yet" << endl; | << "EvtAmp::contract not written yet" << endl; | ||||
▲ Show 20 Lines • Show All 55 Lines • ▼ Show 20 Lines | for ( i = 0; i < allloop[_nontrivial - 1]; i++ ) { | ||||
index++; | index++; | ||||
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) << endl; | EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) << endl; | ||||
} | } | ||||
} | } | ||||
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) | EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) | ||||
<< "-----------------------------------" << endl; | << "-----------------------------------" << endl; | ||||
} | } | ||||
/* | |||||
void EvtAmp::vertex( const EvtComplex& c ) | void EvtAmp::vertex( const EvtComplex& c ) | ||||
{ | { | ||||
int list[1]; | int list[1]; | ||||
list[0] = 0; | list[0] = 0; | ||||
setAmp( list, c ); | setAmp( list, c ); | ||||
} | } | ||||
void EvtAmp::vertex( int i, const EvtComplex& c ) | void EvtAmp::vertex( int i, const EvtComplex& c ) | ||||
Show All 19 Lines | void EvtAmp::vertex( int i, int j, int k, const EvtComplex& c ) | ||||
list[2] = k; | list[2] = k; | ||||
setAmp( list, c ); | setAmp( list, c ); | ||||
} | } | ||||
void EvtAmp::vertex( int* i1, const EvtComplex& c ) | void EvtAmp::vertex( int* i1, const EvtComplex& c ) | ||||
{ | { | ||||
setAmp( i1, c ); | setAmp( i1, c ); | ||||
} | } | ||||
*/ | |||||
EvtAmp& EvtAmp::operator=( const EvtAmp& amp ) | EvtAmp& EvtAmp::operator=( const EvtAmp& amp ) | ||||
{ | { | ||||
int i; | int i; | ||||
_ndaug = amp._ndaug; | _ndaug = amp._ndaug; | ||||
_pstates = amp._pstates; | _pstates = amp._pstates; | ||||
for ( i = 0; i < _ndaug; i++ ) { | for ( i = 0; i < _ndaug; i++ ) { | ||||
Show All 19 Lines |