Line data Source code
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 0 : EvtSpinDensity::EvtSpinDensity(const EvtSpinDensity& density){
34 0 : dim=0;
35 0 : rho=0;
36 :
37 : int i,j;
38 0 : setDim(density.dim);
39 :
40 0 : for(i=0;i<dim;i++){
41 0 : for(j=0;j<dim;j++){
42 0 : rho[i][j]=density.rho[i][j];
43 : }
44 : }
45 0 : }
46 :
47 : EvtSpinDensity& EvtSpinDensity::operator=(const EvtSpinDensity& density){
48 : int i,j;
49 0 : setDim(density.dim);
50 :
51 0 : for(i=0;i<dim;i++){
52 0 : for(j=0;j<dim;j++){
53 0 : rho[i][j]=density.rho[i][j];
54 : }
55 : }
56 :
57 0 : return *this;
58 :
59 : }
60 :
61 0 : EvtSpinDensity::~EvtSpinDensity(){
62 :
63 0 : if (dim!=0){
64 : int i;
65 0 : for(i=0;i<dim;i++) delete [] rho[i];
66 0 : }
67 :
68 0 : delete [] rho;
69 :
70 0 : }
71 :
72 0 : EvtSpinDensity::EvtSpinDensity(){
73 0 : dim=0;
74 0 : rho=0;
75 0 : }
76 :
77 : void EvtSpinDensity::setDim(int n){
78 0 : if (dim==n) return;
79 0 : if (dim!=0){
80 : int i;
81 0 : for(i=0;i<dim;i++) delete [] rho[i];
82 0 : delete [] rho;
83 0 : rho=0;
84 0 : dim=0;
85 0 : }
86 0 : if (n==0) return;
87 0 : dim=n;
88 0 : rho=new EvtComplexPtr[n];
89 : int i;
90 0 : for(i=0;i<n;i++){
91 0 : rho[i]=new EvtComplex[n];
92 : }
93 :
94 :
95 0 : }
96 :
97 : int EvtSpinDensity::getDim() const {
98 0 : return dim;
99 : }
100 :
101 : void EvtSpinDensity::set(int i,int j,const EvtComplex& rhoij){
102 0 : assert(i<dim&&j<dim);
103 0 : rho[i][j]=rhoij;
104 0 : }
105 :
106 : const EvtComplex& EvtSpinDensity::get(int i,int j)const{
107 0 : assert(i<dim&&j<dim);
108 0 : return rho[i][j];
109 : }
110 :
111 : void EvtSpinDensity::setDiag(int n){
112 0 : setDim(n);
113 : int i,j;
114 :
115 0 : for(i=0;i<n;i++){
116 0 : for(j=0;j<n;j++){
117 0 : rho[i][j]=EvtComplex(0.0);
118 : }
119 0 : rho[i][i]=EvtComplex(1.0);
120 : }
121 0 : }
122 :
123 : double EvtSpinDensity::normalizedProb(const EvtSpinDensity& d){
124 :
125 : int i,j;
126 0 : EvtComplex prob(0.0,0.0);
127 : double norm=0.0;
128 :
129 0 : if (dim!=d.dim) {
130 0 : report(Severity::Error,"EvtGen")<<"Not matching dimensions in NormalizedProb"<<endl;
131 0 : ::abort();
132 : }
133 :
134 0 : for(i=0;i<dim;i++){
135 0 : norm+=real(rho[i][i]);
136 0 : for(j=0;j<dim;j++){
137 0 : prob+=rho[i][j]*d.rho[i][j];
138 : }
139 : }
140 :
141 0 : if (imag(prob)>0.00000001*real(prob)) {
142 0 : report(Severity::Error,"EvtGen")<<"Imaginary probability:"<<prob<<" "<<norm<<endl;
143 0 : }
144 0 : if (real(prob)<0.0) {
145 0 : report(Severity::Error,"EvtGen")<<"Negative probability:"<<prob<<" "<<norm<<endl;
146 0 : }
147 :
148 0 : return real(prob)/norm;
149 :
150 0 : }
151 :
152 : int EvtSpinDensity::check(){
153 :
154 0 : if (dim<1) {
155 0 : report(Severity::Error,"EvtGen")<<"dim="<<dim<<"in SpinDensity::Check"<<endl;
156 0 : }
157 :
158 : int i,j;
159 :
160 : double trace(0.0);
161 :
162 0 : for (i=0;i<dim;i++) {
163 0 : trace += abs(rho[i][i]);
164 : }
165 :
166 0 : for(i=0;i<dim;i++){
167 :
168 0 : if (real(rho[i][i])<0.0) return 0;
169 0 : if (imag(rho[i][i])*1000000.0>trace) {
170 0 : report(Severity::Info,"EvtGen") << *this << endl;
171 0 : report(Severity::Info,"EvtGen") << trace << endl;
172 0 : report(Severity::Info,"EvtGen") << "Failing 1"<<endl;
173 0 : return 0;
174 : }
175 : }
176 :
177 0 : for(i=0;i<dim;i++){
178 0 : for(j=i+1;j<dim;j++){
179 0 : if (fabs(real(rho[i][j]-rho[j][i]))>
180 0 : 0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
181 0 : report(Severity::Info,"EvtGen") << "Failing 2"<<endl;
182 0 : return 0;
183 : }
184 0 : if (fabs(imag(rho[i][j]+rho[j][i]))>
185 0 : 0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
186 0 : report(Severity::Info,"EvtGen") << "Failing 3"<<endl;
187 0 : return 0;
188 : }
189 : }
190 : }
191 :
192 0 : return 1;
193 0 : }
194 :
195 : ostream& operator<<(ostream& s,const EvtSpinDensity& d){
196 :
197 : int i,j;
198 :
199 0 : s << endl;
200 0 : s << "Dimension:"<<d.dim<<endl;
201 :
202 0 : for (i=0;i<d.dim;i++){
203 0 : for (j=0;j<d.dim;j++){
204 0 : s << d.rho[i][j]<<" ";
205 : }
206 0 : s <<endl;
207 : }
208 :
209 0 : return s;
210 :
211 : }
212 :
213 :
214 :
|
