]>
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 | 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) | |
48 | ||
49 | }; | |
50 | ||
51 | //Cint craps out here, we protect this part | |
52 | #if !defined(__CINT__) && !defined(__MAKECINT__) | |
53 | ||
54 | //#define DEBUG | |
55 | ||
56 | //______________________________________________________________________________ | |
57 | template <Int_t N> template <class T> inline void AliAODRedCov<N>::GetCovMatrix(T *cmat) const | |
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) { | |
66 | cmat[i*(i+1)/2+j] = (fDiag[j] >= 0. && fDiag[i] >= 0.) ? fODia[(i-1)*i/2+j]*fDiag[j]*fDiag[i]: -999.; | |
67 | #ifdef DEBUG | |
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]); | |
70 | #endif | |
71 | } | |
72 | ||
73 | // Diagonal elements | |
74 | cmat[i*(i+1)/2+i] = (fDiag[i] >= 0.) ? fDiag[i]*fDiag[i] : -999.; | |
75 | #ifdef DEBUG | |
76 | printf("cmat[%2d] = fDiag[%2d]*fDiag[%2d] = %f\n", | |
77 | i*(i+1)/2+i,i,i,cmat[i*(i+1)/2+i]); | |
78 | #endif | |
79 | } | |
80 | } | |
81 | ||
82 | ||
83 | //______________________________________________________________________________ | |
84 | template <Int_t N> template <class T> inline void AliAODRedCov<N>::SetCovMatrix(T *cmat) | |
85 | { | |
86 | // | |
87 | // Sets the external cov matrix | |
88 | // | |
89 | ||
90 | if(cmat) { | |
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 | ||
98 | // Diagonal elements first | |
99 | for(Int_t i=0; i<N; ++i) { | |
100 | fDiag[i] = (cmat[i*(i+1)/2+i] >= 0.) ? TMath::Sqrt(cmat[i*(i+1)/2+i]) : -999.; | |
101 | #ifdef DEBUG | |
102 | printf("fDiag[%2d] = TMath::Sqrt(cmat[%2d]) = %f\n", | |
103 | i,i*(i+1)/2+i, fDiag[i]); | |
104 | #endif | |
105 | } | |
106 | ||
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) { | |
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 | |
120 | #ifdef DEBUG | |
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]); | |
128 | #endif | |
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 | } | |
134 | ||
135 | return; | |
136 | } | |
137 | ||
138 | #undef DEBUG | |
139 | ||
140 | #endif | |
141 | #endif |