correct for omission
[u/mrichter/AliRoot.git] / STEER / AliAODRedCov.h
CommitLineData
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 16template <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 60template <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 87template <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