00001
00005 #include <bpm/bpm_messages.h>
00006 #include <bpm/bpm_analysis.h>
00007 #include <bpm/bpm_nr.h>
00008
00009 int ana_get_svd_coeffs( bpmproc_t **proc, int num_bpms, int num_svd,
00010 int total_num_evts, double *coeffs, int mode ) {
00011
00012
00013
00014
00015 int num_pars = num_bpms;
00016
00017 if (ANA_SVD_TILT & mode) num_pars += num_bpms - 1;
00018
00019
00020
00021 int evt, ibpm, idx, count = 0, good_evt;
00022
00023
00024 gsl_matrix *A = gsl_matrix_calloc( num_svd, num_pars );
00025 gsl_matrix *V = gsl_matrix_calloc(num_pars, num_pars);
00026 gsl_vector *S = gsl_vector_calloc(num_pars);
00027 gsl_vector *work = gsl_vector_calloc(num_pars);
00028 gsl_vector *b = gsl_vector_calloc(num_svd);
00029 gsl_vector *x = gsl_vector_calloc(num_pars);
00030
00031 for ( evt = 0; evt < total_num_evts; evt++ ) {
00032
00033 idx = 0;
00034 good_evt = 1;
00035
00036 for ( ibpm = 1; ibpm < num_bpms; ibpm++ ) {
00037
00038 if ( ana_cutfn( &(proc[ibpm][evt]) ) == BPM_GOOD_EVENT ) {
00039
00040 gsl_matrix_set( A, count, idx++, proc[ibpm][evt].ddc_pos );
00041
00042 if (ANA_SVD_TILT & mode)
00043 gsl_matrix_set( A, count, idx++, proc[ibpm][evt].ddc_slope );
00044 }
00045 else good_evt = 0;
00046 }
00047
00048 if (!good_evt) continue;
00049
00050
00051 gsl_matrix_set( A, count, idx, 1 );
00052 gsl_vector_set( b, count, proc[0][evt].ddc_pos );
00053 count++;
00054 }
00055
00056
00057
00058
00059 gsl_linalg_SV_decomp ( A, V, S, work);
00060 gsl_linalg_SV_solve (A, V, S, b, x);
00061
00062
00063 for (idx = 0; idx < num_pars; idx++) {
00064 coeffs[idx] = gsl_vector_get(x, idx);
00065 }
00066
00067 return BPM_SUCCESS;
00068 }