EvtSpinDensity.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtSpinDensity.cc
//
// Description: Class to reperesent spindensity matrices.
//
// Modification history:
//
//    RYD       May 29,1997       Module created
//
//------------------------------------------------------------------------
//
#include "EvtSpinDensity.hh"
#include "EvtComplex.hh"
#include "EvtPatches.hh"
#include "EvtReport.hh"
#include <assert.h>
#include <iostream>
#include <math.h>
#include <stdlib.h>
using std::endl;
using std::ostream;
 
EvtSpinDensity::EvtSpinDensity( const EvtSpinDensity& density ) {
  dim = 0;
  rho = 0;
 
  int i, j;
  SetDim( density.dim );
 
  for ( i = 0; i < dim; i++ )
  {
    for ( j = 0; j < dim; j++ ) { rho[i][j] = density.rho[i][j]; }
  }
}
 
EvtSpinDensity& EvtSpinDensity::operator=( const EvtSpinDensity& density ) {
  int i, j;
  SetDim( density.dim );
 
  for ( i = 0; i < dim; i++ )
  {
    for ( j = 0; j < dim; j++ ) { rho[i][j] = density.rho[i][j]; }
  }
 
  return *this;
}
 
EvtSpinDensity::~EvtSpinDensity() {
  if ( dim != 0 )
  {
    int i;
    for ( i = 0; i < dim; i++ ) delete[] rho[i];
  }
 
  delete[] rho;
}
 
EvtSpinDensity::EvtSpinDensity() {
  dim = 0;
  rho = 0;
}
 
void EvtSpinDensity::SetDim( int n ) {
  if ( dim == n ) return;
  if ( dim != 0 )
  {
    int i;
    for ( i = 0; i < dim; i++ ) delete[] rho[i];
    delete[] rho;
    rho = 0;
    dim = 0;
  }
  if ( n == 0 ) return;
  dim = n;
  rho = new EvtComplexPtr[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 ) {
  SetDim( n );
  int i, j;
 
  for ( i = 0; i < n; i++ )
  {
    for ( j = 0; j < n; j++ ) { rho[i][j] = EvtComplex( 0.0 ); }
    rho[i][i] = EvtComplex( 1.0 );
  }
}
 
double EvtSpinDensity::NormalizedProb( const EvtSpinDensity& d ) {
 
  int        i, j;
  EvtComplex prob( 0.0, 0.0 );
  double     norm = 0.0;
 
  if ( dim != d.dim )
  {
    report( ERROR, "EvtGen" ) << "Not matching dimensions in NormalizedProb" << endl;
    ::abort();
  }
 
  for ( i = 0; i < dim; i++ )
  {
    norm += real( rho[i][i] );
    for ( j = 0; j < dim; j++ ) { prob += rho[i][j] * d.rho[i][j]; }
  }
 
  if ( imag( prob ) > 0.00000001 * real( prob ) )
  { report( ERROR, "EvtGen" ) << "Imaginary probability:" << prob << " " << norm << endl; }
  if ( real( prob ) < 0.0 )
  { report( ERROR, "EvtGen" ) << "Negative probability:" << prob << " " << norm << endl; }
 
  return real( prob ) / norm;
}
 
int EvtSpinDensity::Check() {
 
  if ( dim < 1 )
  { report( ERROR, "EvtGen" ) << "dim=" << dim << "in SpinDensity::Check" << endl; }
 
  int i, j;
 
  for ( i = 0; i < dim; i++ )
  {
 
    if ( real( rho[i][i] ) < 0.0 ) return 0;
    if ( imag( rho[i][i] ) * 1000000.0 > abs( rho[i][i] ) )
    {
      report( INFO, "EvtGen" ) << "Failing 1" << endl;
      return 0;
    }
  }
 
  for ( i = 0; i < dim; i++ )
  {
    for ( j = i + 1; j < dim; j++ )
    {
      if ( fabs( real( rho[i][j] - rho[j][i] ) ) >
           0.00000001 * ( abs( rho[i][i] ) + abs( rho[j][j] ) ) )
      {
        report( INFO, "EvtGen" ) << "Failing 2" << endl;
        return 0;
      }
      if ( fabs( imag( rho[i][j] + rho[j][i] ) ) >
           0.00000001 * ( abs( rho[i][i] ) + abs( rho[j][j] ) ) )
      {
        report( INFO, "EvtGen" ) << "Failing 3" << endl;
        return 0;
      }
    }
  }
 
  return 1;
}
 
ostream& operator<<( ostream& s, const EvtSpinDensity& d ) {
 
  int i, j;
 
  s << endl;
  s << "Dimension:" << d.dim << endl;
 
  for ( i = 0; i < d.dim; i++ )
  {
    for ( j = 0; j < d.dim; j++ ) { s << d.rho[i][j] << " "; }
    s << endl;
  }
 
  return s;
}