Page MenuHomeHEPForge
Differential D87 Diff 389 src/EvtGenBase/EvtSpinDensity.cpp

Changeset View

Standalone View

View Options

src/EvtGenBase/EvtSpinDensity.cpp

Show All 18 Lines
***********************************************************************/ ***********************************************************************/
#include "EvtGenBase/EvtSpinDensity.hh" #include "EvtGenBase/EvtSpinDensity.hh"
#include "EvtGenBase/EvtComplex.hh" #include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtReport.hh"
#include <assert.h> #include <cstring>
#include <iostream> #include <iostream>
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
using std::endl; using std::endl;
using std::ostream; using std::ostream;
EvtSpinDensity::EvtSpinDensity( const EvtSpinDensity& density ) EvtSpinDensity::EvtSpinDensity( const EvtSpinDensity& density )
{ {
dim = 0; dim = 0;
rho = nullptr; rho = nullptr;
int i, j;
setDim( density.dim ); setDim( density.dim );
// memmove(rho,density.rho,dim*dim*sizeof(EvtComplex));
for ( i = 0; i < dim; i++ ) { for ( int i = 0; i < dim * dim; i++ )
for ( j = 0; j < dim; j++ ) { rho[i] = density.rho[i];
rho[i][j] = density.rho[i][j];
}
}
} }
EvtSpinDensity& EvtSpinDensity::operator=( const EvtSpinDensity& density ) EvtSpinDensity& EvtSpinDensity::operator=( const EvtSpinDensity& density )
{ {
int i, j;
setDim( density.dim ); setDim( density.dim );
// memmove(rho,density.rho,dim*dim*sizeof(EvtComplex));
for ( i = 0; i < dim; i++ ) { for ( int i = 0; i < dim * dim; i++ )
for ( j = 0; j < dim; j++ ) { rho[i] = density.rho[i];
rho[i][j] = density.rho[i][j];
}
}
return *this; return *this;
} }
EvtSpinDensity::~EvtSpinDensity() EvtSpinDensity::~EvtSpinDensity()
{ {
if ( dim != 0 ) { if ( dim != 0 )
int i;
for ( i = 0; i < dim; i++ )
delete[] rho[i];
}
delete[] rho; delete[] rho;
} }
EvtSpinDensity::EvtSpinDensity()
{
dim = 0;
rho = nullptr;
}
void EvtSpinDensity::setDim( int n ) void EvtSpinDensity::setDim( int n )
{ {
if ( dim == n ) if ( dim == n )
return; return;
if ( dim != 0 ) { if ( dim != 0 ) {
int i;
for ( i = 0; i < dim; i++ )
delete[] rho[i];
delete[] rho; delete[] rho;
rho = nullptr; rho = nullptr;
dim = 0; dim = 0;
} }
if ( n == 0 ) if ( n == 0 )
return; return;
dim = n; dim = n;
rho = new EvtComplexPtr[n]; rho = new EvtComplex[n * n];
int i;
for ( i = 0; i < n; i++ ) {
rho[i] = new EvtComplex[n];
}
}
int EvtSpinDensity::getDim() const
{
return dim;
}
void EvtSpinDensity::set( int i, int j, const EvtComplex& rhoij )
{
assert( i < dim && j < dim );
rho[i][j] = rhoij;
}
const EvtComplex& EvtSpinDensity::get( int i, int j ) const
{
assert( i < dim && j < dim );
return rho[i][j];
} }
void EvtSpinDensity::setDiag( int n ) void EvtSpinDensity::setDiag( int n )
{ {
setDim( n ); setDim( n );
int i, j; for ( int i = 0; i < n * n; i++ )
rho[i] = 0;
for ( i = 0; i < n; i++ ) { for ( int i = 0; i < n; i++ )
for ( j = 0; j < n; j++ ) { rho[i * dim + i] = 1.0;
rho[i][j] = EvtComplex( 0.0 );
}
rho[i][i] = EvtComplex( 1.0 );
}
} }
double EvtSpinDensity::normalizedProb( const EvtSpinDensity& d ) double EvtSpinDensity::normalizedProb( const EvtSpinDensity& d ) const
{ {
int i, j;
EvtComplex prob( 0.0, 0.0 );
double norm = 0.0;
if ( dim != d.dim ) { if ( dim != d.dim ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Not matching dimensions in NormalizedProb" << endl; << "Not matching dimensions in NormalizedProb" << endl;
::abort(); ::abort();
} }
for ( i = 0; i < dim; i++ ) { double norm = 0.0;
norm += real( rho[i][i] ); for ( int i = 0; i < dim; i++ )
for ( j = 0; j < dim; j++ ) { norm += real( rho[i * dim + i] );
prob += rho[i][j] * d.rho[i][j];
} EvtComplex prob;
} for ( int i = 0, imax = dim * dim; i < imax; i++ )
prob += rho[i] * d.rho[i];
if ( imag( prob ) > 0.00000001 * real( prob ) ) { if ( imag( prob ) > 0.00000001 * real( prob ) ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Imaginary probability:" << prob << " " << norm << endl; << "Imaginary probability:" << prob << " " << norm << endl;
} }
if ( real( prob ) < 0.0 ) { if ( real( prob ) < 0.0 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Negative probability:" << prob << " " << norm << endl; << "Negative probability:" << prob << " " << norm << endl;
Show All 9 Lines if ( dim < 1 ) {
<< "dim=" << dim << "in SpinDensity::Check" << endl; << "dim=" << dim << "in SpinDensity::Check" << endl;
} }
int i, j; int i, j;
double trace( 0.0 ); double trace( 0.0 );
for ( i = 0; i < dim; i++ ) { for ( i = 0; i < dim; i++ ) {
trace += abs( rho[i][i] ); trace += abs( rho[i * dim + i] );
} }
for ( i = 0; i < dim; i++ ) { for ( i = 0; i < dim; i++ ) {
if ( real( rho[i][i] ) < 0.0 ) if ( real( rho[i * dim + i] ) < 0.0 )
return 0; return 0;
if ( imag( rho[i][i] ) * 1000000.0 > trace ) { if ( imag( rho[i * dim + i] ) * 1000000.0 > trace ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << *this << endl; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << *this << endl;
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << trace << endl; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << trace << endl;
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 1" << endl; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 1" << endl;
return 0; return 0;
} }
} }
for ( i = 0; i < dim; i++ ) { for ( i = 0; i < dim; i++ ) {
for ( j = i + 1; j < dim; j++ ) { for ( j = i + 1; j < dim; j++ ) {
if ( fabs( real( rho[i][j] - rho[j][i] ) ) > if ( fabs( real( rho[i * dim + j] - rho[j * dim + i] ) ) >
0.00000001 * ( abs( rho[i][i] ) + abs( rho[j][j] ) ) ) { 0.00000001 *
( abs( rho[i * dim + i] ) + abs( rho[j * dim + j] ) ) ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 2" << endl; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 2" << endl;
return 0; return 0;
} }
if ( fabs( imag( rho[i][j] + rho[j][i] ) ) > if ( fabs( imag( rho[i * dim + j] + rho[j * dim + i] ) ) >
0.00000001 * ( abs( rho[i][i] ) + abs( rho[j][j] ) ) ) { 0.00000001 *
( abs( rho[i * dim + i] ) + abs( rho[j * dim + j] ) ) ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 3" << endl; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 3" << endl;
return 0; return 0;
} }
} }
} }
return 1; return 1;
} }
ostream& operator<<( ostream& s, const EvtSpinDensity& d ) ostream& operator<<( ostream& s, const EvtSpinDensity& d )
{ {
int i, j;
s << endl; s << endl;
s << "Dimension:" << d.dim << endl; s << "Dimension:" << d.dim << endl;
for ( i = 0; i < d.dim; i++ ) { for ( int i = 0; i < d.dim; i++ ) {
for ( j = 0; j < d.dim; j++ ) { for ( int j = 0; j < d.dim; j++ ) {
s << d.rho[i][j] << " "; s << d.get( i, j ) << " ";
} }
s << endl; s << endl;
} }
return s; return s;
} }