X-Git-Url: http://git.treefish.org/~alex/phys/u1casc.git/blobdiff_plain/a47125a5c94fc75125547d64c693e7384ed9b70d..607713e7bb60eb68ad4285a805f10943a7cec2ff:/u1casc-ordinary/obs_diagcorr.hpp

diff --git a/u1casc-ordinary/obs_diagcorr.hpp b/u1casc-ordinary/obs_diagcorr.hpp
index f51df7c..a8d88b5 100644
--- a/u1casc-ordinary/obs_diagcorr.hpp
+++ b/u1casc-ordinary/obs_diagcorr.hpp
@@ -46,9 +46,7 @@ private:
 
   gsl_vector_complex **measurements_disco;
   
-  static void cdiag (gsl_vector *v, gsl_matrix_complex *m);
-
-  static void invertcm (gsl_matrix_complex *invm, const gsl_matrix_complex *m);
+  static void cdiag (gsl_vector *v, gsl_matrix_complex *evec, gsl_matrix_complex *m);
 };
 
 obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr", 
@@ -97,23 +95,8 @@ void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) {
   }
 };
 
-void obs_diagcorr::invertcm (gsl_matrix_complex *invm, const gsl_matrix_complex *m)
-{
-  int s;
-  gsl_permutation * perm = gsl_permutation_alloc (m->size1);
-  
-  gsl_matrix_complex *mtmp = gsl_matrix_complex_alloc(m->size1,m->size2);
-  gsl_matrix_complex_memcpy(mtmp, m);
-
-  gsl_linalg_complex_LU_decomp (mtmp, perm, &s);
-  gsl_linalg_complex_LU_invert (mtmp, perm, invm);
-
-  gsl_matrix_complex_free(mtmp);
-}
-
-void obs_diagcorr::cdiag (gsl_vector *v, gsl_matrix_complex *m)
+void obs_diagcorr::cdiag (gsl_vector *v, gsl_matrix_complex *evec, gsl_matrix_complex *m)
 {
-  gsl_matrix_complex *evec = gsl_matrix_complex_alloc(4,4);  
   gsl_eigen_hermv_workspace *wspace = gsl_eigen_hermv_alloc(100);
 
   gsl_eigen_hermv(m, v, evec, wspace);
@@ -133,7 +116,8 @@ void obs_diagcorr::_finish() {
   gsl_vector *jackres = gsl_vector_alloc(4);
   gsl_vector *jackerrornorm = gsl_vector_alloc(4);
   gsl_vector_complex *totaldisco = gsl_vector_complex_alloc(4);
-  gsl_matrix_complex *invC0 = NULL;
+  gsl_matrix_complex *evecres = gsl_matrix_complex_alloc(4,4);
+  gsl_matrix_complex *evecerrornorm = gsl_matrix_complex_alloc(4,4);
   
   gsl_vector_complex_set_zero(totaldisco);
   for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
@@ -142,6 +126,7 @@ void obs_diagcorr::_finish() {
   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
     gsl_matrix_complex_set_zero(totalval);
     gsl_vector_set_zero(jackerrornorm);
+    gsl_matrix_complex_set_zero(evecerrornorm);
 
     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
       gsl_matrix_complex_add( totalval, measurements[imeas][itsep] );
@@ -160,15 +145,7 @@ void obs_diagcorr::_finish() {
 						 )
 				);
       }
-
-    if (invC0 == NULL) {
-      invC0 = gsl_matrix_complex_alloc(4,4);
-      invertcm(invC0, tmpmatrix); 
-    }
-
-    gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,0.0),
-		    invC0, tmpmatrix, gsl_complex_rect(0.0,0.0), tmpmatrix2);
-    cdiag(jackres, tmpmatrix2);
+    cdiag(jackres, evecres, tmpmatrix);
 
     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++) {
       gsl_matrix_complex_memcpy (tmpmatrix, totalval);
@@ -187,31 +164,43 @@ void obs_diagcorr::_finish() {
 						   discopart
 						   )
 				  );
-	}
-
-      gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,0.0),
-		      invC0, tmpmatrix, gsl_complex_rect(0.0,0.0), tmpmatrix2);
+	}      
+      cdiag(tmpvec, tmpmatrix2, tmpmatrix);
       
-      cdiag(tmpvec, tmpmatrix2);
       gsl_vector_sub(tmpvec, jackres);
       gsl_vector_memcpy(tmpvec2, tmpvec);
       gsl_vector_mul(tmpvec, tmpvec2);
       gsl_vector_add(jackerrornorm, tmpvec);
+
+      gsl_matrix_complex_sub(tmpmatrix2, evecres);
+      for (int i=0; i<4; i++)
+	for (int j=0; j<4; j++) {
+	  gsl_complex tmperr =
+	    gsl_complex_rect( pow(GSL_REAL( gsl_matrix_complex_get( tmpmatrix2, i, j ) ),2),
+			      pow(GSL_IMAG( gsl_matrix_complex_get( tmpmatrix2, i, j ) ),2));
+	  gsl_matrix_complex_set(tmpmatrix, i, j, tmperr);
+	}
+      gsl_matrix_complex_add( evecerrornorm, tmpmatrix );
     } 
     gsl_vector_scale( jackerrornorm,
 		      (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas );
-
+    gsl_matrix_complex_scale( evecerrornorm, gsl_complex_rect( (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas, 0.0 ) );
+    
     *out << O815->paraQ->getParaVals();
     *out << "\t" << itsep;
     for (int iev = 0; iev < 4; iev++) {
       *out << "\t" << gsl_vector_get(jackres,iev)
 	   << "\t" << sqrt(gsl_vector_get(jackerrornorm,iev));
+      for (int iinter=0; iinter<4; iinter++)
+	*out << "\t" << GSL_REAL(gsl_matrix_complex_get(evecres, iinter, iev))
+	     << "\t" << sqrt( GSL_REAL(gsl_matrix_complex_get(evecerrornorm, iinter, iev)) );
     }
     *out << endl;
   }
 
   gsl_matrix_complex_free(totalval);
   gsl_matrix_complex_free(tmpmatrix);
+  gsl_matrix_complex_free(evecres);
   gsl_vector_free(tmpvec);
   gsl_vector_complex_free(tmpvecc);
   gsl_vector_free(tmpvec2);