]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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 | for(i=0;i<dim;i++){ | |
161 | ||
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; | |
165 | return 0; | |
166 | } | |
167 | } | |
168 | ||
169 | for(i=0;i<dim;i++){ | |
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; | |
174 | return 0; | |
175 | } | |
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; | |
179 | return 0; | |
180 | } | |
181 | } | |
182 | } | |
183 | ||
184 | return 1; | |
185 | } | |
186 | ||
187 | ostream& operator<<(ostream& s,const EvtSpinDensity& d){ | |
188 | ||
189 | int i,j; | |
190 | ||
191 | s << endl; | |
192 | s << "Dimension:"<<d.dim<<endl; | |
193 | ||
194 | for (i=0;i<d.dim;i++){ | |
195 | for (j=0;j<d.dim;j++){ | |
196 | s << d.rho[i][j]<<" "; | |
197 | } | |
198 | s <<endl; | |
199 | } | |
200 | ||
201 | return s; | |
202 | ||
203 | } | |
204 | ||
205 | ||
206 |