Fix for definitions for CINT
[u/mrichter/AliRoot.git] / 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