]>
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: |
38 | AliAODRedCov() {} | |
39 | virtual ~AliAODRedCov() {} | |
40 | template <class T> void GetCovMatrix(T *cmat) const; | |
41 | template <class T> void SetCovMatrix(T *cmat); | |
42 | ||
43 | private: | |
44 | Double32_t fDiag[N]; // Diagonal elements | |
45 | Double32_t fODia[N*(N-1)/2]; // [-1, 1,8] 8 bit precision for off diagonal elements | |
46 | ||
47 | ClassDef(AliAODRedCov,1) | |
df9db588 | 48 | |
4d209fca | 49 | }; |
df9db588 | 50 | |
c98a50a1 | 51 | //Cint craps out here, we protect this part |
4d209fca | 52 | #if !defined(__CINT__) && !defined(__MAKECINT__) |
c98a50a1 | 53 | |
b4116488 | 54 | //#define DEBUG |
c98a50a1 | 55 | |
4d209fca | 56 | //______________________________________________________________________________ |
c98a50a1 | 57 | template <Int_t N> template <class T> inline void AliAODRedCov<N>::GetCovMatrix(T *cmat) const |
4d209fca | 58 | { |
59 | // | |
60 | // Returns the external cov matrix | |
61 | // | |
62 | ||
63 | for(Int_t i=0; i<N; ++i) { | |
64 | // Off diagonal elements | |
65 | for(Int_t j=0; j<i; ++j) { | |
31fd97b2 | 66 | cmat[i*(i+1)/2+j] = (fDiag[j] >= 0. && fDiag[i] >= 0.) ? fODia[(i-1)*i/2+j]*fDiag[j]*fDiag[i]: -999.; |
c98a50a1 | 67 | #ifdef DEBUG |
31fd97b2 | 68 | printf("cmat[%2d] = fODia[%2d]*fDiag[%2d]*fDiag[%2d] = %f\n", |
69 | i*(i+1)/2+j,(i-1)*i/2+j,j,i,cmat[i*(i+1)/2+j]); | |
c98a50a1 | 70 | #endif |
31fd97b2 | 71 | } |
4d209fca | 72 | |
73 | // Diagonal elements | |
31fd97b2 | 74 | cmat[i*(i+1)/2+i] = (fDiag[i] >= 0.) ? fDiag[i]*fDiag[i] : -999.; |
c98a50a1 | 75 | #ifdef DEBUG |
31fd97b2 | 76 | printf("cmat[%2d] = fDiag[%2d]*fDiag[%2d] = %f\n", |
77 | i*(i+1)/2+i,i,i,cmat[i*(i+1)/2+i]); | |
c98a50a1 | 78 | #endif |
4d209fca | 79 | } |
80 | } | |
df9db588 | 81 | |
df9db588 | 82 | |
4d209fca | 83 | //______________________________________________________________________________ |
c98a50a1 | 84 | template <Int_t N> template <class T> inline void AliAODRedCov<N>::SetCovMatrix(T *cmat) |
4d209fca | 85 | { |
86 | // | |
87 | // Sets the external cov matrix | |
88 | // | |
df9db588 | 89 | |
4d209fca | 90 | if(cmat) { |
31fd97b2 | 91 | |
92 | #ifdef DEBUG | |
93 | for (Int_t i=0; i<(N*(N+1))/2; i++) { | |
94 | printf("cmat[%d] = %f\n", i, cmat[i]); | |
95 | } | |
96 | #endif | |
97 | ||
4d209fca | 98 | // Diagonal elements first |
99 | for(Int_t i=0; i<N; ++i) { | |
31fd97b2 | 100 | fDiag[i] = (cmat[i*(i+1)/2+i] >= 0.) ? TMath::Sqrt(cmat[i*(i+1)/2+i]) : -999.; |
c98a50a1 | 101 | #ifdef DEBUG |
31fd97b2 | 102 | printf("fDiag[%2d] = TMath::Sqrt(cmat[%2d]) = %f\n", |
103 | i,i*(i+1)/2+i, fDiag[i]); | |
c98a50a1 | 104 | #endif |
31fd97b2 | 105 | } |
106 | ||
4d209fca | 107 | // ... then the ones off diagonal |
108 | for(Int_t i=0; i<N; ++i) | |
109 | // Off diagonal elements | |
110 | for(Int_t j=0; j<i; ++j) { | |
31fd97b2 | 111 | fODia[(i-1)*i/2+j] = (fDiag[i] > 0. && fDiag[j] > 0.) ? cmat[i*(i+1)/2+j]/(fDiag[j]*fDiag[i]) : 0.; |
112 | // check for division by zero (due to diagonal element of 0) and for fDiag != -999. (due to negative input diagonal element). | |
113 | if (fODia[(i-1)*i/2+j]>1.) { // check upper boundary | |
114 | #ifdef DEBUG | |
115 | printf("out of bounds: %f\n", fODia[(i-1)*i/2+j]); | |
116 | #endif | |
117 | fODia[(i-1)*i/2+j] = 1.; | |
118 | } | |
119 | if (fODia[(i-1)*i/2+j]<-1.) { // check lower boundary | |
c98a50a1 | 120 | #ifdef DEBUG |
31fd97b2 | 121 | printf("out of bounds: %f\n", fODia[(i-1)*i/2+j]); |
122 | #endif | |
123 | fODia[(i-1)*i/2+j] = -1.; | |
124 | } | |
125 | #ifdef DEBUG | |
126 | printf("fODia[%2d] = cmat[%2d]/(fDiag[%2d]*fDiag[%2d]) = %f\n", | |
127 | (i-1)*i/2+j,i*(i+1)/2+j,j,i,fODia[(i-1)*i/2+j]); | |
c98a50a1 | 128 | #endif |
4d209fca | 129 | } |
130 | } else { | |
131 | for(Int_t i=0; i< N; ++i) fDiag[i]=-999.; | |
132 | for(Int_t i=0; i< N*(N-1)/2; ++i) fODia[i]=0.; | |
133 | } | |
31fd97b2 | 134 | |
135 | return; | |
4d209fca | 136 | } |
137 | ||
c98a50a1 | 138 | #undef DEBUG |
139 | ||
4d209fca | 140 | #endif |
df9db588 | 141 | #endif |