TEvtGen/EvtGenBase/EvtSpinDensity.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtSpinDensity.cc
  12 //
  13 // Description: Class to reperesent spindensity matrices.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD       May 29,1997       Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include <iostream>
  24 #include <math.h>
  25 #include <assert.h>
  26 #include "EvtGenBase/EvtComplex.hh"
  27 #include "EvtGenBase/EvtSpinDensity.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 using std::endl;
  30 using std::ostream;
  31
  32
  33 EvtSpinDensity::EvtSpinDensity(const EvtSpinDensity& density){
  34   dim=0;
  35   rho=0;
  36
  37   int i,j;
  38   setDim(density.dim);
  39
  40   for(i=0;i<dim;i++){
  41     for(j=0;j<dim;j++){
  42       rho[i][j]=density.rho[i][j];
  43     }
  44   }
  45 }
  46
  47 EvtSpinDensity& EvtSpinDensity::operator=(const EvtSpinDensity& density){
  48   int i,j;
  49   setDim(density.dim);
  50
  51   for(i=0;i<dim;i++){
  52     for(j=0;j<dim;j++){
  53       rho[i][j]=density.rho[i][j];
  54     }
  55   }
  56
  57   return *this;
  58
  59 }
  60
  61 EvtSpinDensity::~EvtSpinDensity(){
  62
  63   if (dim!=0){
  64     int i;
  65     for(i=0;i<dim;i++) delete [] rho[i];
  66   }
  67
  68   delete [] rho;
  69
  70 }
  71
  72 EvtSpinDensity::EvtSpinDensity(){
  73   dim=0;
  74   rho=0;
  75 }
  76
  77 void EvtSpinDensity::setDim(int n){
  78   if (dim==n) return;
  79   if (dim!=0){
  80     int i;
  81     for(i=0;i<dim;i++) delete [] rho[i];
  82     delete [] rho;
  83     rho=0;
  84     dim=0;
  85   }
  86   if (n==0) return;
  87   dim=n;
  88   rho=new EvtComplexPtr[n];
  89   int i;
  90   for(i=0;i<n;i++){
  91     rho[i]=new EvtComplex[n];
  92   }
  93
  94
  95 }
  96
  97 int EvtSpinDensity::getDim() const {
  98   return dim;
  99 }
 100
 101 void EvtSpinDensity::set(int i,int j,const EvtComplex& rhoij){
 102   assert(i<dim&&j<dim);
 103   rho[i][j]=rhoij;
 104 }
 105
 106 const EvtComplex& EvtSpinDensity::get(int i,int j)const{
 107   assert(i<dim&&j<dim);
 108   return rho[i][j];
 109 }
 110
 111 void EvtSpinDensity::setDiag(int n){
 112   setDim(n);
 113   int i,j;
 114
 115   for(i=0;i<n;i++){
 116     for(j=0;j<n;j++){
 117       rho[i][j]=EvtComplex(0.0);
 118     }
 119     rho[i][i]=EvtComplex(1.0);
 120   }
 121 }
 122
 123 double EvtSpinDensity::normalizedProb(const EvtSpinDensity& d){
 124
 125   int i,j;
 126   EvtComplex prob(0.0,0.0);
 127   double norm=0.0;
 128
 129   if (dim!=d.dim) {
 130     report(ERROR,"EvtGen")<<"Not matching dimensions in NormalizedProb"<<endl;
 131     ::abort();
 132   }
 133
 134   for(i=0;i<dim;i++){
 135     norm+=real(rho[i][i]);
 136     for(j=0;j<dim;j++){
 137       prob+=rho[i][j]*d.rho[i][j];
 138     }
 139   }
 140
 141   if (imag(prob)>0.00000001*real(prob)) {
 142     report(ERROR,"EvtGen")<<"Imaginary probability:"<<prob<<" "<<norm<<endl;
 143   }
 144   if (real(prob)<0.0) {
 145     report(ERROR,"EvtGen")<<"Negative probability:"<<prob<<" "<<norm<<endl;
 146   }
 147
 148   return real(prob)/norm;
 149
 150 }
 151
 152 int EvtSpinDensity::check(){
 153
 154   if (dim<1) {
 155     report(ERROR,"EvtGen")<<"dim="<<dim<<"in SpinDensity::Check"<<endl;
 156   }
 157
 158   int i,j;
 159
 160   double trace(0.0);
 161
 162   for (i=0;i<dim;i++) {
 163     trace += abs(rho[i][i]);
 164   }
 165
 166   for(i=0;i<dim;i++){
 167
 168     if (real(rho[i][i])<0.0) return 0;
 169     if (imag(rho[i][i])*1000000.0>trace) {
 170       report(INFO,"EvtGen") << *this << endl;
 171       report(INFO,"EvtGen") << trace << endl;
 172       report(INFO,"EvtGen") << "Failing 1"<<endl;
 173       return 0;
 174     }
 175   }
 176
 177   for(i=0;i<dim;i++){
 178     for(j=i+1;j<dim;j++){
 179       if (fabs(real(rho[i][j]-rho[j][i]))>
 180           0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
 181         report(INFO,"EvtGen") << "Failing 2"<<endl;
 182         return 0;
 183       }
 184       if (fabs(imag(rho[i][j]+rho[j][i]))>
 185           0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
 186         report(INFO,"EvtGen") << "Failing 3"<<endl;
 187         return 0;
 188       }
 189     }
 190   }
 191
 192   return 1;
 193 }
 194
 195 ostream& operator<<(ostream& s,const EvtSpinDensity& d){
 196
 197   int i,j;
 198
 199   s << endl;
 200   s << "Dimension:"<<d.dim<<endl;
 201
 202   for (i=0;i<d.dim;i++){
 203     for (j=0;j<d.dim;j++){
 204      s << d.rho[i][j]<<" ";
 205     }
 206     s <<endl;
 207   }
 208
 209   return s;
 210
 211 }
 212
 213
 214