1 //--------------------------------------------------------------------------
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.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtSpinDensity.cc
13 // Description: Class to reperesent spindensity matrices.
15 // Modification history:
17 // RYD May 29,1997 Module created
19 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
26 #include "EvtGenBase/EvtComplex.hh"
27 #include "EvtGenBase/EvtSpinDensity.hh"
28 #include "EvtGenBase/EvtReport.hh"
33 EvtSpinDensity::EvtSpinDensity(const EvtSpinDensity& density){
42 rho[i][j]=density.rho[i][j];
47 EvtSpinDensity& EvtSpinDensity::operator=(const EvtSpinDensity& density){
53 rho[i][j]=density.rho[i][j];
61 EvtSpinDensity::~EvtSpinDensity(){
65 for(i=0;i<dim;i++) delete [] rho[i];
72 EvtSpinDensity::EvtSpinDensity(){
77 void EvtSpinDensity::setDim(int n){
81 for(i=0;i<dim;i++) delete [] rho[i];
88 rho=new EvtComplexPtr[n];
91 rho[i]=new EvtComplex[n];
97 int EvtSpinDensity::getDim() const {
101 void EvtSpinDensity::set(int i,int j,const EvtComplex& rhoij){
102 assert(i<dim&&j<dim);
106 const EvtComplex& EvtSpinDensity::get(int i,int j)const{
107 assert(i<dim&&j<dim);
111 void EvtSpinDensity::setDiag(int n){
117 rho[i][j]=EvtComplex(0.0);
119 rho[i][i]=EvtComplex(1.0);
123 double EvtSpinDensity::normalizedProb(const EvtSpinDensity& d){
126 EvtComplex prob(0.0,0.0);
130 report(ERROR,"EvtGen")<<"Not matching dimensions in NormalizedProb"<<endl;
135 norm+=real(rho[i][i]);
137 prob+=rho[i][j]*d.rho[i][j];
141 if (imag(prob)>0.00000001*real(prob)) {
142 report(ERROR,"EvtGen")<<"Imaginary probability:"<<prob<<" "<<norm<<endl;
144 if (real(prob)<0.0) {
145 report(ERROR,"EvtGen")<<"Negative probability:"<<prob<<" "<<norm<<endl;
148 return real(prob)/norm;
152 int EvtSpinDensity::check(){
155 report(ERROR,"EvtGen")<<"dim="<<dim<<"in SpinDensity::Check"<<endl;
162 if (real(rho[i][i])<0.0) return 0;
163 if (imag(rho[i][i])*1000000.0>abs(rho[i][i])) {
164 report(INFO,"EvtGen") << "Failing 1"<<endl;
170 for(j=i+1;j<dim;j++){
171 if (fabs(real(rho[i][j]-rho[j][i]))>
172 0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
173 report(INFO,"EvtGen") << "Failing 2"<<endl;
176 if (fabs(imag(rho[i][j]+rho[j][i]))>
177 0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
178 report(INFO,"EvtGen") << "Failing 3"<<endl;
187 ostream& operator<<(ostream& s,const EvtSpinDensity& d){
192 s << "Dimension:"<<d.dim<<endl;
194 for (i=0;i<d.dim;i++){
195 for (j=0;j<d.dim;j++){
196 s << d.rho[i][j]<<" ";