Line data Source code
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 0 : 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 660 : AliAODRedCov() {
39 2058 : for(Int_t i=0; i<N; i++) fDiag[i] = 0.;
40 4272 : for(Int_t i=0; i<N*(N-1)/2; i++) fODia[i] = 0.;
41 330 : }
42 990 : 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 832 : 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 72 : for(Int_t i=0; i<N; ++i) {
67 : // Off diagonal elements
68 96 : for(Int_t j=0; j<i; ++j) {
69 96 : 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 72 : 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 8 : }
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 322 : 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 2014 : for(Int_t i=0; i<N; ++i) {
103 2538 : 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 2014 : for(Int_t i=0; i<N; ++i)
112 : // Off diagonal elements
113 5562 : for(Int_t j=0; j<i; ++j) {
114 7644 : 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 1935 : 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 4 : fODia[(i-1)*i/2+j] = 1.;
121 4 : }
122 1935 : 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 17 : fODia[(i-1)*i/2+j] = -1.;
127 17 : }
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 161 : } else {
134 0 : for(Int_t i=0; i< N; ++i) fDiag[i]=-999.;
135 0 : for(Int_t i=0; i< N*(N-1)/2; ++i) fODia[i]=0.;
136 : }
137 :
138 161 : return;
139 : }
140 :
141 : #undef DEBUG
142 :
143 : #endif
144 : #endif
|
