3 /* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * See cxx source for full Copyright notice */
8 //-------------------------------------------------------------------------
11 //-------------------------------------------------------------------------
16 template <Int_t N> class AliAODRedCov {
20 // Class containing reduced cov matrix, see example here for a track
26 // Y fOdia[ 0] fDiag[ 1]
28 // Z fOdia[ 1] fOdia[ 2] fDiag[ 2]
30 // Px fOdia[ 3] fOdia[ 4] fOdia[ 5] fDiag[ 3]
32 // Py fOdia[ 6] fOdia[ 7] fOdia[ 8] fOdia[ 9] fDiag[ 4]
34 // Pz fOdia[10] fOdia[11] fOdia[12] fOdia[13] fOdia[14] fDiag[ 5]
39 virtual ~AliAODRedCov() {}
40 template <class T> void GetCovMatrix(T *cmat) const;
41 template <class T> void SetCovMatrix(T *cmat);
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
47 ClassDef(AliAODRedCov,1)
51 //Cint craps out here, we protect this part
52 #if !defined(__CINT__) && !defined(__MAKECINT__)
56 //______________________________________________________________________________
57 template <Int_t N> template <class T> inline void AliAODRedCov<N>::GetCovMatrix(T *cmat) const
60 // Returns the external cov matrix
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.;
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]);
74 cmat[i*(i+1)/2+i] = (fDiag[i] >= 0.) ? fDiag[i]*fDiag[i] : -999.;
76 printf("cmat[%2d] = fDiag[%2d]*fDiag[%2d] = %f\n",
77 i*(i+1)/2+i,i,i,cmat[i*(i+1)/2+i]);
83 //______________________________________________________________________________
84 template <Int_t N> template <class T> inline void AliAODRedCov<N>::SetCovMatrix(T *cmat)
87 // Sets the external cov matrix
93 for (Int_t i=0; i<(N*(N+1))/2; i++) {
94 printf("cmat[%d] = %f\n", i, cmat[i]);
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.;
102 printf("fDiag[%2d] = TMath::Sqrt(cmat[%2d]) = %f\n",
103 i,i*(i+1)/2+i, fDiag[i]);
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
115 printf("out of bounds: %f\n", fODia[(i-1)*i/2+j]);
117 fODia[(i-1)*i/2+j] = 1.;
119 if (fODia[(i-1)*i/2+j]<-1.) { // check lower boundary
121 printf("out of bounds: %f\n", fODia[(i-1)*i/2+j]);
123 fODia[(i-1)*i/2+j] = -1.;
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]);
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.;