#ifndef AliAODRedCov_H
#define AliAODRedCov_H
#include <Rtypes.h>
#include <TMath.h>
template <Int_t N> class AliAODRedCov {
public:
AliAODRedCov() {
for(Int_t i=0; i<N; i++) fDiag[i] = 0.;
for(Int_t i=0; i<N*(N-1)/2; i++) fODia[i] = 0.;
}
virtual ~AliAODRedCov() {}
template <class T> void GetCovMatrix(T *cmat) const;
template <class T> void SetCovMatrix(T *cmat);
private:
Double32_t fDiag[N];
Double32_t fODia[N*(N-1)/2];
ClassDef(AliAODRedCov,1)
};
#if !defined(__CINT__) && !defined(__MAKECINT__)
template <Int_t N> template <class T> inline void AliAODRedCov<N>::GetCovMatrix(T *cmat) const
{
for(Int_t i=0; i<N; ++i) {
for(Int_t j=0; j<i; ++j) {
cmat[i*(i+1)/2+j] = (fDiag[j] >= 0. && fDiag[i] >= 0.) ? fODia[(i-1)*i/2+j]*fDiag[j]*fDiag[i]: -999.;
#ifdef DEBUG
printf("cmat[%2d] = fODia[%2d]*fDiag[%2d]*fDiag[%2d] = %f\n",
i*(i+1)/2+j,(i-1)*i/2+j,j,i,cmat[i*(i+1)/2+j]);
#endif
}
cmat[i*(i+1)/2+i] = (fDiag[i] >= 0.) ? fDiag[i]*fDiag[i] : -999.;
#ifdef DEBUG
printf("cmat[%2d] = fDiag[%2d]*fDiag[%2d] = %f\n",
i*(i+1)/2+i,i,i,cmat[i*(i+1)/2+i]);
#endif
}
}
template <Int_t N> template <class T> inline void AliAODRedCov<N>::SetCovMatrix(T *cmat)
{
if(cmat) {
#ifdef DEBUG
for (Int_t i=0; i<(N*(N+1))/2; i++) {
printf("cmat[%d] = %f\n", i, cmat[i]);
}
#endif
for(Int_t i=0; i<N; ++i) {
fDiag[i] = (cmat[i*(i+1)/2+i] >= 0.) ? TMath::Sqrt(cmat[i*(i+1)/2+i]) : -999.;
#ifdef DEBUG
printf("fDiag[%2d] = TMath::Sqrt(cmat[%2d]) = %f\n",
i,i*(i+1)/2+i, fDiag[i]);
#endif
}
for(Int_t i=0; i<N; ++i)
for(Int_t j=0; j<i; ++j) {
fODia[(i-1)*i/2+j] = (fDiag[i] > 0. && fDiag[j] > 0.) ? cmat[i*(i+1)/2+j]/(fDiag[j]*fDiag[i]) : 0.;
if (fODia[(i-1)*i/2+j]>1.) {
#ifdef DEBUG
printf("out of bounds: %f\n", fODia[(i-1)*i/2+j]);
#endif
fODia[(i-1)*i/2+j] = 1.;
}
if (fODia[(i-1)*i/2+j]<-1.) {
#ifdef DEBUG
printf("out of bounds: %f\n", fODia[(i-1)*i/2+j]);
#endif
fODia[(i-1)*i/2+j] = -1.;
}
#ifdef DEBUG
printf("fODia[%2d] = cmat[%2d]/(fDiag[%2d]*fDiag[%2d]) = %f\n",
(i-1)*i/2+j,i*(i+1)/2+j,j,i,fODia[(i-1)*i/2+j]);
#endif
}
} else {
for(Int_t i=0; i< N; ++i) fDiag[i]=-999.;
for(Int_t i=0; i< N*(N-1)/2; ++i) fODia[i]=0.;
}
return;
}
#undef DEBUG
#endif
#endif