STEER/AOD/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