]>
Commit | Line | Data |
---|---|---|
1 | #ifndef AliAODRedCov_H | |
2 | #define AliAODRedCov_H | |
3 | /* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. * | |
4 | * See cxx source for full Copyright notice */ | |
5 | ||
6 | /* $Id$ */ | |
7 | ||
8 | //------------------------------------------------------------------------- | |
9 | // Reduced Cov Matrix | |
10 | // Author: fca | |
11 | //------------------------------------------------------------------------- | |
12 | ||
13 | #include <Rtypes.h> | |
14 | #include <TMath.h> | |
15 | ||
16 | template <Int_t N> class AliAODRedCov { | |
17 | ||
18 | ||
19 | // | |
20 | // Class containing reduced cov matrix, see example here for a track | |
21 | // | |
22 | // X Y Z Px Py Pz | |
23 | // | |
24 | // X fDiag[ 0] | |
25 | // | |
26 | // Y fOdia[ 0] fDiag[ 1] | |
27 | // | |
28 | // Z fOdia[ 1] fOdia[ 2] fDiag[ 2] | |
29 | // | |
30 | // Px fOdia[ 3] fOdia[ 4] fOdia[ 5] fDiag[ 3] | |
31 | // | |
32 | // Py fOdia[ 6] fOdia[ 7] fOdia[ 8] fOdia[ 9] fDiag[ 4] | |
33 | // | |
34 | // Pz fOdia[10] fOdia[11] fOdia[12] fOdia[13] fOdia[14] fDiag[ 5] | |
35 | // | |
36 | ||
37 | public: | |
38 | AliAODRedCov() { | |
39 | for(Int_t i=0; i<N; i++) fDiag[i] = 0.; | |
40 | for(Int_t i=0; i<N*(N-1)/2; i++) fODia[i] = 0.; | |
41 | } | |
42 | virtual ~AliAODRedCov() {} | |
43 | template <class T> void GetCovMatrix(T *cmat) const; | |
44 | template <class T> void SetCovMatrix(T *cmat); | |
45 | ||
46 | private: | |
47 | Double32_t fDiag[N]; // Diagonal elements | |
48 | Double32_t fODia[N*(N-1)/2]; // [-1, 1,8] 8 bit precision for off diagonal elements | |
49 | ||
50 | ClassDef(AliAODRedCov,1) | |
51 | ||
52 | }; | |
53 | ||
54 | //Cint craps out here, we protect this part | |
55 | #if !defined(__CINT__) && !defined(__MAKECINT__) | |
56 | ||
57 | //#define DEBUG | |
58 | ||
59 | //______________________________________________________________________________ | |
60 | template <Int_t N> template <class T> inline void AliAODRedCov<N>::GetCovMatrix(T *cmat) const | |
61 | { | |
62 | // | |
63 | // Returns the external cov matrix | |
64 | // | |
65 | ||
66 | for(Int_t i=0; i<N; ++i) { | |
67 | // Off diagonal elements | |
68 | for(Int_t j=0; j<i; ++j) { | |
69 | cmat[i*(i+1)/2+j] = (fDiag[j] >= 0. && fDiag[i] >= 0.) ? fODia[(i-1)*i/2+j]*fDiag[j]*fDiag[i]: -999.; | |
70 | #ifdef DEBUG | |
71 | printf("cmat[%2d] = fODia[%2d]*fDiag[%2d]*fDiag[%2d] = %f\n", | |
72 | i*(i+1)/2+j,(i-1)*i/2+j,j,i,cmat[i*(i+1)/2+j]); | |
73 | #endif | |
74 | } | |
75 | ||
76 | // Diagonal elements | |
77 | cmat[i*(i+1)/2+i] = (fDiag[i] >= 0.) ? fDiag[i]*fDiag[i] : -999.; | |
78 | #ifdef DEBUG | |
79 | printf("cmat[%2d] = fDiag[%2d]*fDiag[%2d] = %f\n", | |
80 | i*(i+1)/2+i,i,i,cmat[i*(i+1)/2+i]); | |
81 | #endif | |
82 | } | |
83 | } | |
84 | ||
85 | ||
86 | //______________________________________________________________________________ | |
87 | template <Int_t N> template <class T> inline void AliAODRedCov<N>::SetCovMatrix(T *cmat) | |
88 | { | |
89 | // | |
90 | // Sets the external cov matrix | |
91 | // | |
92 | ||
93 | if(cmat) { | |
94 | ||
95 | #ifdef DEBUG | |
96 | for (Int_t i=0; i<(N*(N+1))/2; i++) { | |
97 | printf("cmat[%d] = %f\n", i, cmat[i]); | |
98 | } | |
99 | #endif | |
100 | ||
101 | // Diagonal elements first | |
102 | for(Int_t i=0; i<N; ++i) { | |
103 | fDiag[i] = (cmat[i*(i+1)/2+i] >= 0.) ? TMath::Sqrt(cmat[i*(i+1)/2+i]) : -999.; | |
104 | #ifdef DEBUG | |
105 | printf("fDiag[%2d] = TMath::Sqrt(cmat[%2d]) = %f\n", | |
106 | i,i*(i+1)/2+i, fDiag[i]); | |
107 | #endif | |
108 | } | |
109 | ||
110 | // ... then the ones off diagonal | |
111 | for(Int_t i=0; i<N; ++i) | |
112 | // Off diagonal elements | |
113 | for(Int_t j=0; j<i; ++j) { | |
114 | fODia[(i-1)*i/2+j] = (fDiag[i] > 0. && fDiag[j] > 0.) ? cmat[i*(i+1)/2+j]/(fDiag[j]*fDiag[i]) : 0.; | |
115 | // check for division by zero (due to diagonal element of 0) and for fDiag != -999. (due to negative input diagonal element). | |
116 | if (fODia[(i-1)*i/2+j]>1.) { // check upper boundary | |
117 | #ifdef DEBUG | |
118 | printf("out of bounds: %f\n", fODia[(i-1)*i/2+j]); | |
119 | #endif | |
120 | fODia[(i-1)*i/2+j] = 1.; | |
121 | } | |
122 | if (fODia[(i-1)*i/2+j]<-1.) { // check lower boundary | |
123 | #ifdef DEBUG | |
124 | printf("out of bounds: %f\n", fODia[(i-1)*i/2+j]); | |
125 | #endif | |
126 | fODia[(i-1)*i/2+j] = -1.; | |
127 | } | |
128 | #ifdef DEBUG | |
129 | printf("fODia[%2d] = cmat[%2d]/(fDiag[%2d]*fDiag[%2d]) = %f\n", | |
130 | (i-1)*i/2+j,i*(i+1)/2+j,j,i,fODia[(i-1)*i/2+j]); | |
131 | #endif | |
132 | } | |
133 | } else { | |
134 | for(Int_t i=0; i< N; ++i) fDiag[i]=-999.; | |
135 | for(Int_t i=0; i< N*(N-1)/2; ++i) fODia[i]=0.; | |
136 | } | |
137 | ||
138 | return; | |
139 | } | |
140 | ||
141 | #undef DEBUG | |
142 | ||
143 | #endif | |
144 | #endif |