00001
00005 #include <bpm/bpm_messages.h>
00006 #include <bpm/bpm_analysis.h>
00007
00008 int ana_compute_residual( bpmproc_t **proc, int num_bpms, int num_evts,
00009 double *coeffs, int mode, double *mean, double *rms ) {
00010
00011
00012 double *good_res = (double*) calloc( num_evts, sizeof(double));
00013 int evt, ibpm, idx, count = 0, good_evt;
00014 double residual;
00015
00016 for ( evt = 0; evt < num_evts; evt++ ) {
00017
00018 idx = 0;
00019 good_evt = 1;
00020 residual = proc[0][evt].ddc_pos;
00021
00022 for ( ibpm = 1; ibpm < num_bpms; ibpm++ ) {
00023
00024 if ( ana_cutfn( &(proc[ibpm][evt]) ) == BPM_GOOD_EVENT ) {
00025
00026 residual -= proc[ibpm][evt].ddc_pos * coeffs[idx++];
00027
00028 if (ANA_SVD_TILT & mode)
00029 residual -= proc[ibpm][evt].ddc_slope * coeffs[idx++];
00030 }
00031 else good_evt = 0;
00032 }
00033
00034 if (!good_evt) continue;
00035
00036
00037 good_res[count] = residual - coeffs[idx];
00038 count++;
00039 }
00040
00041
00042 *mean = 0;
00043 for (evt = 0; evt < count; evt++) {
00044 *mean += good_res[evt];
00045 }
00046
00047 *mean /= (double) count;
00048
00049
00050 *rms = 0;
00051 for (evt = 0; evt < count; evt++) {
00052
00053 *rms += (good_res[evt] - *mean) * (good_res[evt] - *mean);
00054 }
00055
00056 *rms = sqrt( *rms / count );
00057
00058 free(good_res);
00059
00060 return BPM_SUCCESS;
00061 }