]>
Commit | Line | Data |
---|---|---|
0e318744 | 1 | #ifndef AliAODRedCov_H |
2 | #define AliAODRedCov_H | |
df9db588 | 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> | |
4d209fca | 14 | #include <TMath.h> |
df9db588 | 15 | |
0e318744 | 16 | template <Int_t N> class AliAODRedCov { |
df9db588 | 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 | ||
4d209fca | 37 | public: |
7ba6f91a | 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 | } | |
4d209fca | 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) | |
df9db588 | 51 | |
4d209fca | 52 | }; |
df9db588 | 53 | |
c98a50a1 | 54 | //Cint craps out here, we protect this part |
4d209fca | 55 | #if !defined(__CINT__) && !defined(__MAKECINT__) |
c98a50a1 | 56 | |
b4116488 | 57 | //#define DEBUG |
c98a50a1 | 58 | |
4d209fca | 59 | //______________________________________________________________________________ |
c98a50a1 | 60 | template <Int_t N> template <class T> inline void AliAODRedCov<N>::GetCovMatrix(T *cmat) const |
4d209fca | 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) { | |
31fd97b2 | 69 | cmat[i*(i+1)/2+j] = (fDiag[j] >= 0. && fDiag[i] >= 0.) ? fODia[(i-1)*i/2+j]*fDiag[j]*fDiag[i]: -999.; |
c98a50a1 | 70 | #ifdef DEBUG |
31fd97b2 | 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]); | |
c98a50a1 | 73 | #endif |
31fd97b2 | 74 | } |
4d209fca | 75 | |
76 | // Diagonal elements | |
31fd97b2 | 77 | cmat[i*(i+1)/2+i] = (fDiag[i] >= 0.) ? fDiag[i]*fDiag[i] : -999.; |
c98a50a1 | 78 | #ifdef DEBUG |
31fd97b2 | 79 | printf("cmat[%2d] = fDiag[%2d]*fDiag[%2d] = %f\n", |
80 | i*(i+1)/2+i,i,i,cmat[i*(i+1)/2+i]); | |
c98a50a1 | 81 | #endif |
4d209fca | 82 | } |
83 | } | |
df9db588 | 84 | |
df9db588 | 85 | |
4d209fca | 86 | //______________________________________________________________________________ |
c98a50a1 | 87 | template <Int_t N> template <class T> inline void AliAODRedCov<N>::SetCovMatrix(T *cmat) |
4d209fca | 88 | { |
89 | // | |
90 | // Sets the external cov matrix | |
91 | // | |
df9db588 | 92 | |
4d209fca | 93 | if(cmat) { |
31fd97b2 | 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 | ||
4d209fca | 101 | // Diagonal elements first |
102 | for(Int_t i=0; i<N; ++i) { | |
31fd97b2 | 103 | fDiag[i] = (cmat[i*(i+1)/2+i] >= 0.) ? TMath::Sqrt(cmat[i*(i+1)/2+i]) : -999.; |
c98a50a1 | 104 | #ifdef DEBUG |
31fd97b2 | 105 | printf("fDiag[%2d] = TMath::Sqrt(cmat[%2d]) = %f\n", |
106 | i,i*(i+1)/2+i, fDiag[i]); | |
c98a50a1 | 107 | #endif |
31fd97b2 | 108 | } |
109 | ||
4d209fca | 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) { | |
31fd97b2 | 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 | |
c98a50a1 | 123 | #ifdef DEBUG |
31fd97b2 | 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]); | |
c98a50a1 | 131 | #endif |
4d209fca | 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 | } | |
31fd97b2 | 137 | |
138 | return; | |
4d209fca | 139 | } |
140 | ||
c98a50a1 | 141 | #undef DEBUG |
142 | ||
4d209fca | 143 | #endif |
df9db588 | 144 | #endif |