void GetCovarianceMatrix( const TMatrixD &data, TMatrixD &covMat ){
   
   Int_t n = data.GetNrows();
   Int_t d = data.GetNcols();
   Double_t means[5] = {0.0};

   for(Int_t i = 0; i < n; ++i){
      for( Int_t j =0; j < n; ++j){
         means[j] += data[i][j]/(Double_t)n;
      }
   }

   covMat.ResizeTo( d, d ); 
    for( Int_t j =0; j < d; ++j){
        for( Int_t k=0; k < d; k++){
		covMat[j][k] = 0.0;
        }
    }
   
   for(Int_t i = 0; i < n; ++i){
      for( Int_t j =0; j < d; ++j){
         for( Int_t k=0; k < d; k++){
			covMat[j][k] += (means[j] - data[i][j])*(means[k] - data[i][k])/(n - 1.0);
         }
      }
   }         


   delete []means;
}
         

void imatrix() {

   TFile *f = new TFile("tree.root");
   TTree *T = (TTree*)f->Get("dimu");
   const Int_t nl = 5;
   char *leaves[nl] = {"InvMass", "mMet", "m1metDphi", "m1m2Dphi" , "MuCentrali"};
   TLeaf *leaf[nl];
   for (Int_t i=0;i<nl;i++) leaf[i] = T->GetLeaf(leaves[i]);
   
   Int_t N = (Int_t)T->GetEntries();
   TMatrixD matrix(N,nl);
   for (Int_t entry=0;entry<N;entry++) {
      for (Int_t i=0;i<nl;i++) {
         leaf[i]->GetBranch()->GetEntry(entry);
         matrix[entry][i] = leaf[i]->GetValue();
      }
   }
   matrix.Print("Data Matrix");

   TMatrixD covMatrix;

   GetCovarianceMatrix(matrix, covMatrix);
   covMatrix.Print("covMatrix");

   TMatrixD matInv(covMatrix);
   matInv.Invert().Print("matInv");

   //TDecompSVD svd(covMatrix);
   //svd.Decompose();
   //TMatrixD svdInv = svd.Invert();
   //svdInv.Print("svdInv");

   TMatrixD unit = covMatrix*matInv;
   unit.Print("covMatrix*matInv");
}