Number of sigma pedestal cut increased to 4
[u/mrichter/AliRoot.git] / STEER / AliAODRedCov.h
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