TMatrixD get_corr_matr(const TMatrixD matr_cov) {

  TMatrixD matr = matr_cov;
  matr.Zero();
  Int_t n_rows = matr.GetNrows();
  for (Int_t ii = 0; ii < n_rows; ii++) {
    for (Int_t jj = 0; jj < n_rows; jj++) {
      matr(ii,jj) = matr_cov(ii,jj) 
								/ ( sqrt(matr_cov(ii,ii)) * sqrt(matr_cov(jj,jj)) );
    }
  }
  return matr;

}

// ===========================================
void test_2(){


		double fitpar[3] = {3.036010,
																						2.905944,
																						5.582753};

		double errpar[3] = {0.376452,
																						0.157458,
																						1.068356};


		double cov[9] = {0.141846,
																			-0.011840,
																			-0.144900,
																			-0.011840,
																			0.024803,
																			0.158485,
																			-0.144900,
																			0.158485,
																			1.143070};

		int npar=3;

  matr_cov     = TMatrixD(npar,npar,cov);
		par_orig     = TVectorD(npar,fitpar);
		par_orig_err = TVectorD(npar,errpar);
	
		TMatrixD matr_corr    = get_corr_matr(matr_cov); 

		int nPar = npar;

		// calculate the elements of the upper-triangular matrix L 
		// that gives Lt*L = C
    
		// where Lt is the transpose of L (the "square-root method")
  
		L = TMatrixD(nPar,nPar);

		for(Int_t iPar= 0; iPar < nPar; iPar++) {
				// calculate the diagonal term first
				L(iPar,iPar)= matr_corr(iPar,iPar); //covariance(iPar,iPar);

				for(Int_t k= 0; k < iPar; k++) {
						Double_t tmp= L(k,iPar);
						L(iPar,iPar)-= tmp*tmp;
   }

				L(iPar,iPar)= sqrt(L(iPar,iPar));
   // then the off-diagonal terms
   
				for(Int_t jPar= iPar+1; jPar < nPar; jPar++) {
						L(iPar,jPar)= matr_corr(iPar,jPar); //covariance(iPar,jPar);
						for(Int_t k= 0; k < iPar; k++) {
								L(iPar,jPar)-= L(k,iPar)*L(k,jPar);	
						}
						L(iPar,jPar)/= L(iPar,iPar);
				}
		}

		// remember Lt
		//_Lt=  TMatrixD(TMatrixD::kTransposed,L);

		TMatrixD _Lt(TMatrixD::kTransposed,L);

		//L.Print();
		//_Lt.Print();
				


		// create a vector of unit Gaussian variables
  TVectorD g(nPar);
		
		int nrand = 10000;
	

		double *mean = new double[npar];
		double *rms  = new double[npar];
		double *newpar = new double[npar*nrand];

		double par_val;
		
  TRandom coef;
	
		//g.Print();

		//nrand=1;

		/* initilize mean and rms */
		for(int j=0;j<npar;j++){mean[j]=0; rms[j]=0;}
																		
		/* this is the MonteCarlo Loop */
		for(int i=0;i<nrand;i++){

				for(Int_t k= 0; k < nPar; k++) g(k)= coef.Gaus(0,1.0);
				// multiply this vector by Lt to introduce the 
				//appropriate correlations
				//g*= (*_Lt);

				g = _Lt * g;
				
				/*store the values */

				for(int j=0;j<npar;j++){
						
						par_val = fitpar[j]+g(j);

						mean[j]+=par_val/nrand;

						newpar[npar*i+j] = par_val;
				}

		}/* end of Monte Carlo loop */

		/* check correlation factors */
		TGraph *g1 = new TGraph(nrand);
		TGraph *g2 = new TGraph(nrand);
		TGraph *g3 = new TGraph(nrand);
		
		/* compute the rms */
		for(int i=0;i<nrand;i++){
				g1->SetPoint(i,newpar[npar*i+0], newpar[npar*i+1]);
				g2->SetPoint(i,newpar[npar*i+1], newpar[npar*i+2]);
				g3->SetPoint(i,newpar[npar*i+0], newpar[npar*i+2]);
				
				for(int j=0;j<npar;j++)
						rms[j]+=(newpar[npar*i+j] -mean[j]) * (newpar[npar*i+j]-mean[j])/nrand;
		}
		
		TH1F * h =new TH1F("a","a",30,-3,10);
		for (int j=0;j<nrand;j++)
				h->Fill(newpar[npar*j+0]);
		h->Draw();


		/* print out debuggin lines */
		for(int j=0;j<npar;j++)
				cout<<"Par "<<j<<" Mean "<<mean[j]<<" RMS="<<sqrt(rms[j])<<endl;
				

		/* print Correlation factors */
		cout<<"Correlation 0_1 "<<g1->GetCorrelationFactor()<<
				" while it should be "<<matr_corr(0,1)<<endl;

		cout<<"Correlation 1_2 "<<g2->GetCorrelationFactor()<<
				" while it should be "<<matr_corr(1,2)<<endl;

		cout<<"Correlation 0_2 "<<g3->GetCorrelationFactor()<<
				" while it should be "<<matr_corr(0,2)<<endl;

}