seems to be working.

[phys/u1casc.git] / u1casc-ordinary / obs_diagcorr.hpp

diff --git a/u1casc-ordinary/obs_diagcorr.hpp b/u1casc-ordinary/obs_diagcorr.hpp
index 43d2132905f155414594707a08f183e142ef287d..d52a6cdbf57f7b8c42e0b185e1bf3e992b1ed031 100644 (file)
--- a/u1casc-ordinary/obs_diagcorr.hpp
+++ b/u1casc-ordinary/obs_diagcorr.hpp
@@ -16,6 +16,9 @@
 #include <gsl/gsl_complex_math.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_eigen.h>
 #include <gsl/gsl_complex_math.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_eigen.h>

+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_cblas.h>
+#include <gsl/gsl_blas.h>

 
 using namespace std;
 
 
 using namespace std;
 


@@ -43,7 +46,7 @@ private:
 
   gsl_vector_complex **measurements_disco;
   
 
   gsl_vector_complex **measurements_disco;
   

-  static void cdiag (gsl_vector *v, 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", 
 };
 
 obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr", 


@@ -82,37 +85,39 @@ void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) {
   for (int icorr=0; icorr<4; icorr++)
     for (int jcorr=0; jcorr<4; jcorr++)
       for (int itsep=0; itsep<O815->comargs.lsize[1]/2; itsep++) {
   for (int icorr=0; icorr<4; icorr++)
     for (int jcorr=0; jcorr<4; jcorr++)
       for (int itsep=0; itsep<O815->comargs.lsize[1]/2; itsep++) {

-       gsl_complex tmpc;
-       tmpc.dat[0] = OM[icorr*4+jcorr][itsep].real();
-       tmpc.dat[1] = OM[icorr*4+jcorr][itsep].imag();
-       gsl_matrix_complex_set(measurements[nthmeas][itsep], icorr, jcorr, tmpc);
+       gsl_matrix_complex_set(measurements[nthmeas][itsep], icorr, jcorr,
+                              gsl_complex_rect( OM[icorr*4+jcorr][itsep].real(), OM[icorr*4+jcorr][itsep].imag() ));

       }
       }

-
+  

   for (int icorr=0; icorr<4; icorr++) {
   for (int icorr=0; icorr<4; icorr++) {

-    gsl_complex tmpc;
-    tmpc.dat[0] = OM[16][icorr].real();
-    tmpc.dat[1] = OM[16][icorr].imag();
-    gsl_vector_complex_set(measurements_disco[nthmeas], icorr, tmpc);
+    gsl_vector_complex_set(measurements_disco[nthmeas], icorr,
+                          gsl_complex_rect( OM[16][icorr].real(), OM[16][icorr].imag() ) );

   }
 };
 
   }
 };
 

-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_eigen_herm_workspace *wspace = gsl_eigen_herm_alloc(100);
-  gsl_eigen_herm(m, v, wspace);
-  gsl_eigen_herm_free(wspace);
+  gsl_eigen_hermv_workspace *wspace = gsl_eigen_hermv_alloc(100);
+
+  gsl_eigen_hermv(m, v, evec, wspace);
+  gsl_eigen_hermv_sort(v, evec, GSL_EIGEN_SORT_VAL_ASC);
+  
+  gsl_eigen_hermv_free(wspace);

 }
 
 void obs_diagcorr::_finish() {
   gsl_matrix_complex *totalval = gsl_matrix_complex_alloc(4,4);
   gsl_matrix_complex *tmpmatrix = gsl_matrix_complex_alloc(4,4);
 }
 
 void obs_diagcorr::_finish() {
   gsl_matrix_complex *totalval = gsl_matrix_complex_alloc(4,4);
   gsl_matrix_complex *tmpmatrix = gsl_matrix_complex_alloc(4,4);

+  gsl_matrix_complex *tmpmatrix2 = gsl_matrix_complex_alloc(4,4);

   gsl_vector *tmpvec = gsl_vector_alloc(4);
   gsl_vector_complex *tmpvecc = gsl_vector_complex_alloc(4);
   gsl_vector *tmpvec2 = gsl_vector_alloc(4);
   gsl_vector *jackres = gsl_vector_alloc(4);
   gsl_vector *jackerrornorm = gsl_vector_alloc(4);
   gsl_vector_complex *totaldisco = gsl_vector_complex_alloc(4);
   gsl_vector *tmpvec = gsl_vector_alloc(4);
   gsl_vector_complex *tmpvecc = gsl_vector_complex_alloc(4);
   gsl_vector *tmpvec2 = gsl_vector_alloc(4);
   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 *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++)
     gsl_vector_complex_add( totaldisco, measurements_disco[imeas] );
   gsl_vector_complex_set_zero(totaldisco);
   for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
     gsl_vector_complex_add( totaldisco, measurements_disco[imeas] );


@@ -120,17 +125,13 @@ 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);
   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] );
 
     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
       gsl_matrix_complex_add( totalval, measurements[imeas][itsep] );

-      
-    // compute mean //
-    gsl_complex tmpc2;
-    tmpc2.dat[0] = 1.0/O815->comargs.nmeas;
-    tmpc2.dat[1] = 0;
+    

     gsl_matrix_complex_memcpy (tmpmatrix, totalval);
     gsl_matrix_complex_memcpy (tmpmatrix, totalval);

-    gsl_matrix_complex_scale (tmpmatrix, tmpc2);
-
+    gsl_matrix_complex_scale (tmpmatrix, gsl_complex_rect(1.0/O815->comargs.nmeas, 0.0) );

     for (int icorr=0; icorr<4; icorr++)
       for (int jcorr=0; jcorr<4; jcorr++) {
        gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(totaldisco, icorr),
     for (int icorr=0; icorr<4; icorr++)
       for (int jcorr=0; jcorr<4; jcorr++) {
        gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(totaldisco, icorr),


@@ -143,23 +144,14 @@ void obs_diagcorr::_finish() {
                                                 )
                                );
       }
                                                 )
                                );
       }

+    cdiag(jackres, evecres, tmpmatrix);

     
     

-    cdiag(jackres, tmpmatrix);
-    // END
-

     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++) {
     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++) {

-      gsl_complex tmpc;
-      tmpc.dat[0] = 1.0/(O815->comargs.nmeas-1);
-      tmpc.dat[1] = 0;
-      

       gsl_matrix_complex_memcpy (tmpmatrix, totalval);
       gsl_matrix_complex_sub (tmpmatrix, measurements[imeas][itsep]);
       gsl_matrix_complex_memcpy (tmpmatrix, totalval);
       gsl_matrix_complex_sub (tmpmatrix, measurements[imeas][itsep]);

-
-      gsl_matrix_complex_scale ( tmpmatrix, tmpc );
-
+      gsl_matrix_complex_scale ( tmpmatrix, gsl_complex_rect(1.0/(O815->comargs.nmeas-1), 0.0) );

       gsl_vector_complex_memcpy( tmpvecc, totaldisco );
       gsl_vector_complex_sub( tmpvecc, measurements_disco[imeas] );
       gsl_vector_complex_memcpy( tmpvecc, totaldisco );
       gsl_vector_complex_sub( tmpvecc, measurements_disco[imeas] );

-      

       for (int icorr=0; icorr<4; icorr++)
        for (int jcorr=0; jcorr<4; jcorr++) {
          gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(tmpvecc, icorr),
       for (int icorr=0; icorr<4; icorr++)
        for (int jcorr=0; jcorr<4; jcorr++) {
          gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(tmpvecc, icorr),


@@ -172,29 +164,48 @@ void obs_diagcorr::_finish() {
                                                   )
                                  );
        }
                                                   )
                                  );
        }

+      cdiag(tmpvec, tmpmatrix2, tmpmatrix);

       
       

-      cdiag(tmpvec, tmpmatrix);
-

       gsl_vector_sub(tmpvec, jackres);
       gsl_vector_memcpy(tmpvec2, tmpvec);
       gsl_vector_mul(tmpvec, tmpvec2);
       gsl_vector_sub(tmpvec, jackres);
       gsl_vector_memcpy(tmpvec2, tmpvec);
       gsl_vector_mul(tmpvec, tmpvec2);

-

       gsl_vector_add(jackerrornorm, tmpvec);
       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_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));
     *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);
     }
     *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);
+  gsl_vector_free(jackres);
+  gsl_vector_free(jackerrornorm);
+  gsl_vector_complex_free(totaldisco);

 };
 
 void obs_diagcorr::corrCompute()
 };
 
 void obs_diagcorr::corrCompute()