00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <stdlib.h>
00023 #include <iostream>
00024 #include <math.h>
00025 #include <assert.h>
00026 #include "EvtGenBase/EvtComplex.hh"
00027 #include "EvtGenBase/EvtSpinDensity.hh"
00028 #include "EvtGenBase/EvtReport.hh"
00029 using std::endl;
00030 using std::ostream;
00031
00032
00033 EvtSpinDensity::EvtSpinDensity(const EvtSpinDensity& density){
00034 dim=0;
00035 rho=0;
00036
00037 int i,j;
00038 SetDim(density.dim);
00039
00040 for(i=0;i<dim;i++){
00041 for(j=0;j<dim;j++){
00042 rho[i][j]=density.rho[i][j];
00043 }
00044 }
00045 }
00046
00047 EvtSpinDensity& EvtSpinDensity::operator=(const EvtSpinDensity& density){
00048 int i,j;
00049 SetDim(density.dim);
00050
00051 for(i=0;i<dim;i++){
00052 for(j=0;j<dim;j++){
00053 rho[i][j]=density.rho[i][j];
00054 }
00055 }
00056
00057 return *this;
00058
00059 }
00060
00061 EvtSpinDensity::~EvtSpinDensity(){
00062 if (dim!=0){
00063 int i;
00064 for(i=0;i<dim;i++) delete [] rho[i];
00065 }
00066
00067 delete [] rho;
00068 }
00069
00070 EvtSpinDensity::EvtSpinDensity(){
00071 dim=0;
00072 rho=0;
00073 }
00074
00075 void EvtSpinDensity::SetDim(int n){
00076 if (dim==n) return;
00077 if (dim!=0){
00078 int i;
00079 for(i=0;i<dim;i++) delete [] rho[i];
00080 delete [] rho;
00081 rho=0;
00082 dim=0;
00083 }
00084 if (n==0) return;
00085 dim=n;
00086 rho=new EvtComplexPtr[n];
00087 int i;
00088 for(i=0;i<n;i++){
00089 rho[i]=new EvtComplex[n];
00090 }
00091
00092
00093 }
00094
00095 int EvtSpinDensity::GetDim() const {
00096 return dim;
00097 }
00098
00099 void EvtSpinDensity::Set(int i,int j,const EvtComplex& rhoij){
00100 assert(i<dim&&j<dim);
00101 rho[i][j]=rhoij;
00102 }
00103
00104 const EvtComplex& EvtSpinDensity::Get(int i,int j)const{
00105 assert(i<dim&&j<dim);
00106 return rho[i][j];
00107 }
00108
00109 void EvtSpinDensity::SetDiag(int n){
00110 SetDim(n);
00111 int i,j;
00112
00113 for(i=0;i<n;i++){
00114 for(j=0;j<n;j++){
00115 rho[i][j]=EvtComplex(0.0);
00116 }
00117 rho[i][i]=EvtComplex(1.0);
00118 }
00119 }
00120
00121 double EvtSpinDensity::NormalizedProb(const EvtSpinDensity& d){
00122
00123 int i,j;
00124 EvtComplex prob(0.0,0.0);
00125 double norm=0.0;
00126
00127 if (dim!=d.dim) {
00128 report(ERROR,"EvtGen")<<"Not matching dimensions in NormalizedProb"<<endl;
00129 ::abort();
00130 }
00131
00132 for(i=0;i<dim;i++){
00133 norm+=real(rho[i][i]);
00134 for(j=0;j<dim;j++){
00135 prob+=rho[i][j]*d.rho[i][j];
00136 }
00137 }
00138
00139 if (imag(prob)>0.00000001*real(prob)) {
00140 report(ERROR,"EvtGen")<<"Imaginary probability:"<<prob<<" "<<norm<<endl;
00141 }
00142 if (real(prob)<0.0) {
00143 report(ERROR,"EvtGen")<<"Negative probability:"<<prob<<" "<<norm<<endl;
00144 }
00145
00146 return real(prob)/norm;
00147
00148 }
00149
00150 int EvtSpinDensity::Check(){
00151
00152 if (dim<1) {
00153 report(ERROR,"EvtGen")<<"dim="<<dim<<"in SpinDensity::Check"<<endl;
00154 }
00155
00156 int i,j;
00157
00158 for(i=0;i<dim;i++){
00159
00160 if (real(rho[i][i])<0.0) return 0;
00161 if (imag(rho[i][i])*1000000.0>abs(rho[i][i])) {
00162 report(INFO,"EvtGen") << "Failing 1"<<endl;
00163 return 0;
00164 }
00165 }
00166
00167 for(i=0;i<dim;i++){
00168 for(j=i+1;j<dim;j++){
00169 if (fabs(real(rho[i][j]-rho[j][i]))>
00170 0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
00171 report(INFO,"EvtGen") << "Failing 2"<<endl;
00172 return 0;
00173 }
00174 if (fabs(imag(rho[i][j]+rho[j][i]))>
00175 0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
00176 report(INFO,"EvtGen") << "Failing 3"<<endl;
00177 return 0;
00178 }
00179 }
00180 }
00181
00182 return 1;
00183 }
00184
00185 ostream& operator<<(ostream& s,const EvtSpinDensity& d){
00186
00187 int i,j;
00188
00189 s << endl;
00190 s << "Dimension:"<<d.dim<<endl;
00191
00192 for (i=0;i<d.dim;i++){
00193 for (j=0;j<d.dim;j++){
00194 s << d.rho[i][j]<<" ";
00195 }
00196 s <<endl;
00197 }
00198
00199 return s;
00200
00201 }
00202
00203
00204